lconv_intl_refcount

#ifndef NL_DGP_H_
#define NL_DGP_H_

//#include <boost/numeric/conversion/cast.hpp>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_statistics.h>
#include <asserts.h>

using namespace std;

class NL_Dgp {
	public:
		NL_Dgp () {   }; //default constructor
		~NL_Dgp () {   };//default destructor

        //generate TWO threshold AR processes of the second order. INPUT: A 5x1 vector of coefficients for X (alpha_X), a 5x1 vector of coefficients for Y (alpha_Y),
        //and a template function to generate two random errors for both the processes (gen_RAN). OUTPUT: Tx1 matrices (X and Y).
        template <void gen_RAN (double &, double &, const double, const double, const int, unsigned long)>
        static void gen_TAR (Matrix &X, Matrix &Y, const Matrix alpha_X, const Matrix alpha_Y, const double delta, const double rho, const int choose_alt,
		                     unsigned long seed);
		//generate a univariate threshold AR process and a bivariate threshold AR process of the second orders.
        //INPUT: A 5x1 vector of coefficients for X (alpha_X), a 5x2 matrix of coefficients for Y1 and Y2 (alpha_Y),
        //and a template function to generate two random errors for both the processes (gen_RAN). OUTPUT: a Tx1 vector (X) and Tx1 vectors (Y1 and Y2).
        template <void gen_RAN (double &, double &, double &, const Matrix &, const double, const double, const int, unsigned long)>
        static void gen_TAR (Matrix &X, Matrix &Y1, Matrix &Y2, const Matrix alpha_X, const Matrix alpha_Y, const Matrix &delta, const double alpha5,
                             const double rho1, const int choose_alt, unsigned long seed);
        //generate TWO bilinear processes of the second order. INPUT: A 6x1 vector of coefficients for X (alpha_X), a 6x1 vector of coefficients for Y (alpha_Y),
        //and a template function to generate two random errors for both the processes (gen_RAN). OUTPUT: Tx1 matrices (X and Y).
        template <void gen_RAN (double &, double &, const double, const double, const int, unsigned long)>
        static void gen_Bilinear (Matrix &X, Matrix &Y, const Matrix alpha_X, const Matrix alpha_Y, const double delta, const double rho, const int choose_alt,
		                          unsigned long seed);
		//generate a univariate bilinear process and a bivariate bilinear process of the second orders that may have some dependency.
        //INPUT: A 6x1 vector of coefficients for X (alpha_X), a 6x2 matrix of coefficients for Y1 and Y2 (alpha_Y), and a template function
        //to generate random errors for each individual process (gen_RAN). OUTPUT: a Tx1 vector (X) and Tx1 vectors (Y1 and Y2).
        template <void gen_RAN (double &, double &, double &, const Matrix &, const double, const double, const int, unsigned long)>
        static void gen_Bilinear (Matrix &X, Matrix &Y1, Matrix &Y2, const Matrix alpha_X, const Matrix alpha_Y, const Matrix &delta, const double alpha5,
		                          const double rho1, const int choose_alt, unsigned long seed);
        //generate centered skew-normal random error terms. INPUT: delta in (-1,1), a correlation between xi_1_eps and xi_1_eta (rho), an alternative:
        //choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
        //dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: two random error terms (epsilon and eta).
        static void gen_SN (double &epsilon, double &eta, const double delta, const double rho, const int choose_alt, unsigned long seed);
        //generate a trivariate skew-normal random variable. INPUT: a 3x1 vector (delta) in (-1,1)^3, correlations (rho12 and rho13), an alternative:
        //set choose_alt = 0 to generate independent random variables, choose_alt = 1 to generate *epsilon uncorrelated and dependent with *eta1 and *eta2,
        //choose_alt = 2 to generate dependent random variables. OUTPUT: three random variables (epsilon, eta1, and eta2).
        static void gen_TriSN (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double rho12, const double rho13, const int choose_alt,
		                       unsigned long seed);
        //generate two mixtures of N(0,1) random variables. INPUT: a correlation between xi_1_eps and xi_1_eta (rho) and an alternative:
		//choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
        //dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: two random error terms (epsilon and eta).
        static void gen_MN (double &epsilon, double &eta, const double delta, const double rho, const int choose_alt, unsigned long seed);
        //generate three mixtures of N(0,1) random variables. INPUT: a correlation between xi_1_eps and xi_1_eta1 (rho) and an alternative:
        //choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
        //dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: three random error terms (epsilon, eta1, and eta2).
        static void gen_TriMN (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double alpha5, const double rho1, const int choose_alt,
                               unsigned long seed);
        //generate two centered chi-squared random variables. INPUT: a correlation between xi_1_eps and xi_1_eta (rho) and an alternative:
        //choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
        //dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: two random error terms (epsilon and eta).
        static void gen_ChiSq (double &epsilon, double &eta, const double delta, const double rho, const int choose_alt, unsigned long seed);
        //generate three centered chi-squared random variables. INPUT: a correlation between xi_1_eps and xi_1_eta1 (rho) and an alternative:
        //choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
        //dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: three random error terms (epsilon, eta1, and eta2).
        static void gen_TriChiSq (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double alpha5, const double rho1, const int choose_alt,
                                  unsigned long seed);
        //generate two centered Beta random variables. INPUT: a shape parameter (alpha5) -- alpha5 = 0.0001, 0.1, 1., 2. OUTPUT: two random error terms (epsilon and eta).
        static void gen_Beta (double &epsilon, double &eta, const double alpha5, const double rho, const int choose_alt, unsigned long seed);
        //generate three centered Beta random variables. INPUT: a shape parameter (alpha5) -- alpha5 = 0.0001, 0.1, 1., 2.
        //OUTPUT: three random error terms (epsilon, eta1, and eta2).
        static void gen_TriBeta (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double alpha5, const double rho1, const int choose_alt,
                                 unsigned long seed);
        //generate two vectors of random errors (epsilon and eta). INPUT: constants (delta and rho in (-1,1)), an alternative (choose_alt = 0, 1, 2), and a
        //random generator template (gen_RAN: gen_SN, gen_MN, gen_ChiSq or gen_Beta). OUTPUT: 2 vectors
        template <void gen_RAN (double &, double &, const double, const double, const int, unsigned long)>
        static void gen_RANV (Matrix &epsilon, Matrix &eta, const double delta, const double rho, const int choose_alt, unsigned long seed);
        //generate three vectors of random errors (epsilon, eta1 and eta2). INPUT: a vector (delta in (-1,1)^3) and constants (rho12 and rho13 in (-1,1)),
        //an alternative (choose_alt = 0, 1, 2), and a random generator template (gen_RAN: gen_TriSN, gen_TriMN, gen_TriChiSq or gen_TriBeta). OUTPUT: 3 vectors
        template <void gen_RAN (double &, double &, double &, const Matrix &, const double, const double, const int, unsigned long)>
        static void gen_RANV (Matrix &epsilon, Matrix &eta1, Matrix &eta2, const Matrix &delta, const double rho12, const double rho13, const int choose_alt,
                              unsigned long seed);

		//generate a 3x1 normal random vector with a mean zero and a variance-covariance matrix (x).
        static Matrix multi_norm (gsl_matrix *x, unsigned long seed);
        //run multivariate OLS. INPUT: a Tx1 vector of data on the dependent (Y) and a TxN matrix of data on the independent (X).
        //OUTPUT: a Tx1 vector of residuals (resid) and a Nx1 vector of the OLS estimates (slope)
        static void gen_Resid (Matrix &resid, Matrix &slope, const Matrix X, const Matrix Y);
        //estimate the bilinear regression model by the OLS. INPUT: a Tx1 vector of data on X.
        //OUTPUT: a (T-2)x1 vector of residuals (resid) and a 6x1 vector of the OLS estimates (slope)
        static void est_BL (Matrix &resid, Matrix &slope, const Matrix X);
        //estimate the TAR regression model by the OLS. INPUT: a Tx1 vector of data on X.
        //OUTPUT: a (T-2)x1 vector of residuals (resid) and a 5x1 vector of the OLS estimates (slope)
        static void est_TAR (Matrix &resid, Matrix &slope, const Matrix X);
};

//generate a 3x1 normal random vector with a mean zero and a variance-covariance matrix (x).
Matrix NL_Dgp::multi_norm (gsl_matrix *x, unsigned long seed) { //x is a var-cov matrix
    gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    gsl_vector *mean = gsl_vector_calloc (3); //vector of zero means
    gsl_matrix * L = gsl_matrix_calloc (3, 3);
    gsl_matrix_memcpy (L, x); //copy x into L
    gsl_linalg_cholesky_decomp1 (L); //get a Cholesky decomposition matrix
    gsl_vector *xi_vec = gsl_vector_calloc (3);
	gsl_ran_multivariate_gaussian (r, mean, L, xi_vec); //call the multivariate normal random generator
	Matrix result(3, 1);
	result(1) = gsl_vector_get (xi_vec, 0);
	result(2) = gsl_vector_get (xi_vec, 1);
	result(3) = gsl_vector_get (xi_vec, 2);
	gsl_vector_free (mean); //free memory
	gsl_matrix_free (L);
    gsl_vector_free (xi_vec);
    gsl_rng_free (r);
    return result;
}

//generate TWO bilinear processes of the second order. INPUT: A 6x1 vector of coefficients for X (alpha_X), a 6x1 vector of coefficients for Y (alpha_Y),
//and a template function to generate two random errors for both the processes (gen_RAN). OUTPUT: Tx1 matrices (X and Y).
template <void gen_RAN (double &, double &, const double, const double, const int, unsigned long)>
void NL_Dgp::gen_Bilinear (Matrix &X, Matrix &Y, const Matrix alpha_X, const Matrix alpha_Y, const double delta, const double rho,
                           const int choose_alt, unsigned long seed) {
    gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    auto T = X.nRow();
    ASSERT(T == Y.nRow());
    auto B = 200; //burning the first 200 observations
    Matrix X_tmp(T+B, 1), Y_tmp(T+B, 1);
    unsigned long rseed = gsl_rng_get (r); //a random seed
	gen_RAN (X_tmp(1), Y_tmp(1), delta, rho, choose_alt, seed); //generate initial values
	gen_RAN (X_tmp(2), Y_tmp(2), delta, rho, choose_alt, rseed);
	auto epsilon = 0., eta = 0.;
	for (auto t = 3; t <= T+B; ++t) {
		rseed = gsl_rng_get (r); //a random seed
		gen_RAN (epsilon, eta, delta, rho, choose_alt, rseed); //generate random error terms (with different seeds)
		X_tmp(t) = alpha_X(1) + alpha_X(2)*X_tmp(t-1) + alpha_X(3)*pow(X_tmp(t-1),2.) + alpha_X(4)*X_tmp(t-2) + alpha_X(5)*pow(X_tmp(t-2),2.)
		       + alpha_X(6)*X_tmp(t-1)*X_tmp(t-2) + epsilon;
		Y_tmp(t) = alpha_Y(1) + alpha_Y(2)*Y_tmp(t-1) + alpha_Y(3)*pow(Y_tmp(t-1),2.) + alpha_Y(4)*Y_tmp(t-2) + alpha_Y(5)*pow(Y_tmp(t-2),2.)
		       + alpha_Y(6)*Y_tmp(t-1)*Y_tmp(t-2) + eta;
		if (t > B) {
	        X(t-B) = X_tmp(t);
	    	Y(t-B) = Y_tmp(t);
		}
	}
	gsl_rng_free (r); //free memory
}

//generate a univariate bilinear process and a bivariate bilinear process of the second orders that may have some dependency.
//INPUT: A 6x1 vector of coefficients for X (alpha_X), a 6x2 matrix of coefficients for Y1 and Y2 (alpha_Y), and a template function
//to generate random errors for each individual process (gen_RAN). OUTPUT: a Tx1 vector (X) and Tx1 vectors (Y1 and Y2).
template <void gen_RAN (double &, double &, double &, const Matrix &, const double, const double, const int, unsigned long)>
void NL_Dgp::gen_Bilinear (Matrix &X, Matrix &Y1, Matrix &Y2, const Matrix alpha_X, const Matrix alpha_Y, const Matrix &delta,
                           const double alpha5, const double rho1, const int choose_alt, unsigned long seed) {
    gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    auto T = X.nRow();
    ASSERT(T == Y1.nRow() && T == Y2.nRow());
    auto B = 200; //burning the first 200 observations
    Matrix X_tmp(T+B,1), Y_tmp(T+B,2);
    unsigned long rseed = gsl_rng_get (r); //a random seed
    gen_RAN (X_tmp(1), Y_tmp(1,1), Y_tmp(1,2), delta, alpha5, rho1, choose_alt, seed); //generate initial values
    gen_RAN (X_tmp(2), Y_tmp(2,1), Y_tmp(2,2), delta, alpha5, rho1, choose_alt, rseed);
    auto epsilon = 0., eta1 = 0., eta2 = 0.;
    for (auto t = 3; t <= T + B; ++t) {
    	rseed = gsl_rng_get (r); //a random seed
    	gen_RAN (epsilon, eta1, eta2, delta, alpha5, rho1, choose_alt, rseed); //generate random errors (with different seeds)
    	X_tmp(t) = alpha_X(1) + alpha_X(2)*X_tmp(t-1) + alpha_X(3)*pow(X_tmp(t-1),2.) + alpha_X(4)*X_tmp(t-2) + alpha_X(5)*pow(X_tmp(t-2),2.)
    	           + alpha_X(6)*X_tmp(t-1)*X_tmp(t-2) + epsilon;
    	Y_tmp(t,1) = alpha_Y(1,1) + alpha_Y(2,1)*Y_tmp(t-1,1) + alpha_Y(3,1)*pow(Y_tmp(t-1,1),2.) + alpha_Y(4,1)*Y_tmp(t-2,1)
		           + alpha_Y(5,1)*pow(Y_tmp(t-2,1),2.) + alpha_Y(6,1)*Y_tmp(t-1,1)*Y_tmp(t-2,1) + eta1;
		Y_tmp(t,2) = alpha_Y(1,2) + alpha_Y(2,2)*Y_tmp(t-1,2) + alpha_Y(3,2)*pow(Y_tmp(t-1,2),2.) + alpha_Y(4,2)*Y_tmp(t-2,2)
		           + alpha_Y(5,2)*pow(Y_tmp(t-2,2),2.) + alpha_Y(6,2)*Y_tmp(t-1,2)*Y_tmp(t-2,2) + eta2;
		if (t > B) {
	    	X(t-B) = X_tmp(t);
	    	Y1(t-B) = Y_tmp(t,1);
	    	Y2(t-B) = Y_tmp(t,2);
		}
    }
    gsl_rng_free (r); //free memory
}

//generate TWO threshold AR processes of the second order. INPUT: A 5x1 vector of coefficients for X (alpha_X), a 5x1 vector of coefficients for Y (alpha_Y),
//and a template function to generate two random errors for both the processes (gen_RAN). OUTPUT: Tx1 matrices (X and Y).
template <void gen_RAN (double &, double &, const double, const double, const int, unsigned long)>
void NL_Dgp::gen_TAR (Matrix &X, Matrix &Y, const Matrix alpha_X, const Matrix alpha_Y, const double delta, const double rho, const int choose_alt,
                      unsigned long seed) {
    gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    auto T = X.nRow();
    ASSERT(T == Y.nRow());
    auto B = 200; //burning the first 200 observations
    Matrix X_tmp(T+B,1), Y_tmp(T+B,1);
    unsigned long rseed = gsl_rng_get (r); //a random seed
	gen_RAN (X_tmp(1), Y_tmp(1), delta, rho, choose_alt, seed); //generate initial values
	gen_RAN (X_tmp(2), Y_tmp(2), delta, rho, choose_alt, rseed);
	auto epsilon = 0., eta = 0.;
	auto pos = [](double x) { //get the positive part of a real number (x)
    	if (x > 0.)
    		return x;
    	else
    		return 0.;
	};
	for (auto t = 3; t <= T+B; ++t) {
		rseed = gsl_rng_get (r); //a random seed
		gen_RAN (epsilon, eta, delta, rho, choose_alt, rseed); //generate random error terms (with different seeds)
		X_tmp(t) = alpha_X(1) + alpha_X(2)*X_tmp(t-1) + alpha_X(3)*pos(X_tmp(t-1)) + alpha_X(4)*X_tmp(t-2) + alpha_X(5)*pos(X_tmp(t-2)) + epsilon;
		Y_tmp(t) = alpha_Y(1) + alpha_Y(2)*Y_tmp(t-1) + alpha_Y(3)*pos(Y_tmp(t-1)) + alpha_Y(4)*Y_tmp(t-2) + alpha_Y(5)*pos(Y_tmp(t-2)) + eta;
		if (t > B) {
	        X(t-B) = X_tmp(t);
	    	Y(t-B) = Y_tmp(t);
		}
	}
	gsl_rng_free (r); //free memory
}

//generate a univariate threshold AR process and a bivariate threshold AR process of the second orders.
//INPUT: A 5x1 vector of coefficients for X (alpha_X), a 5x2 matrix of coefficients for Y1 and Y2 (alpha_Y),
//and a template function to generate two random errors for both the processes (gen_RAN). OUTPUT: a Tx1 vector (X) and Tx1 vectors (Y1 and Y2).
template <void gen_RAN (double &, double &, double &, const Matrix &, const double, const double, const int, unsigned long)>
void NL_Dgp::gen_TAR (Matrix &X, Matrix &Y1, Matrix &Y2, const Matrix alpha_X, const Matrix alpha_Y, const Matrix &delta, const double alpha5,
                      const double rho1, const int choose_alt, unsigned long seed) {
    gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
	auto T = X.nRow();
    ASSERT(T == Y1.nRow() && T == Y2.nRow());
    auto B = 200; //burning the first 200 observations
    Matrix X_tmp(T+B, 1), Y_tmp(T+B, 2);
    unsigned long rseed = gsl_rng_get (r); //a random seed
    gen_RAN (X_tmp(1), Y_tmp(1,1), Y_tmp(1,2), delta, alpha5, rho1, choose_alt, seed); //generate initial values
    gen_RAN (X_tmp(2), Y_tmp(2,1), Y_tmp(2,2), delta, alpha5, rho1, choose_alt, rseed);
    double epsilon = 0., eta1 = 0., eta2 = 0.;
    auto pos = [](double x) { //get the positive part of a real number (x)
    	if (x > 0.)
    		return x;
    	else
    		return 0.;
	};
    for (auto t = 3; t <= T + B; ++t) {
    	rseed = gsl_rng_get (r); //a random seed
    	gen_RAN (epsilon, eta1, eta2, delta, alpha5, rho1, choose_alt, rseed); //generate random errors (with different seeds)
    	X_tmp(t) = alpha_X(1) + alpha_X(2)*X_tmp(t-1) + alpha_X(3)*pos(X_tmp(t-1)) + alpha_X(4)*X_tmp(t-2) + alpha_X(5)*pos(X_tmp(t-2)) + epsilon;
    	Y_tmp(t,1) = alpha_Y(1,1) + alpha_Y(2,1)*Y_tmp(t-1,1) + alpha_Y(3,1)*pos(Y_tmp(t-1,1)) + alpha_Y(4,1)*Y_tmp(t-2,1) + alpha_Y(5,1)*pos(Y_tmp(t-2,1))
		            + eta1;
		Y_tmp(t,2) = alpha_Y(1,2) + alpha_Y(2,2)*Y_tmp(t-1,2) + alpha_Y(3,2)*pos(Y_tmp(t-1,2)) + alpha_Y(4,2)*Y_tmp(t-2,2) + alpha_Y(5,2)*pos(Y_tmp(t-2,2))
		            + eta2;
		if (t > B) {
	        X(t-B) = X_tmp(t);
	    	Y1(t-B) = Y_tmp(t,1);
	    	Y2(t-B) = Y_tmp(t,2);
		}
    }
	gsl_rng_free (r); //free memory
}

//generate centered skew-normal random error terms. INPUT: delta in (-1,1), a correlation between xi_1_eps and xi_1_eta (rho), an alternative:
//choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
//dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: two random error terms (epsilon and eta).
void NL_Dgp::gen_SN (double &epsilon, double &eta, const double delta, const double rho, const int choose_alt, unsigned long seed) {
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    double xi_0_eps = gsl_ran_ugaussian (r);
    double xi_0_eta = gsl_ran_ugaussian (r);
    double xi_1_eps = gsl_ran_ugaussian (r);
	double xi_1_eta = 0.;
	switch (choose_alt) {
        case 0:
                xi_1_eta = gsl_ran_ugaussian (r); //independent error terms
                break;
        case 1:
        	    gsl_ran_bivariate_gaussian (r, 1., 1., rho, &xi_1_eps, &xi_1_eta); //correlated error terms
        	    break;
        case 2:
        	    if (fabs(xi_1_eps) <= 1.54) //set threshold equal to 1.54 for zero correlation
	                xi_1_eta = xi_1_eps;
	            else
	                xi_1_eta = -xi_1_eps;
	            break;
	    case 3:
	            if (fabs(xi_1_eps) <= 2.) //set threshold equal to 2.0 for non-zero correlation
	                xi_1_eta = xi_1_eps;
	            else
	                xi_1_eta = -xi_1_eps;
	            break;
	    default:
    		    cerr << "NL_Dgp::gen_SM: This choice is not in the switch list. Make sure that your choice is valid!\n";
                exit(0);
	}
	epsilon = delta * fabs(xi_0_eps) + sqrt(1 - pow(delta,2.)) * xi_1_eps - sqrt((double) 2/M_PI) * delta;
	eta =     delta * fabs(xi_0_eta) + sqrt(1 - pow(delta,2.)) * xi_1_eta - sqrt((double) 2/M_PI) * delta;
	gsl_rng_free (r); //free memory
}

//generate a trivariate skew-normal random variable. INPUT: a 3x1 vector (delta) in (-1,1)^3, correlations (rho12 and rho13), an alternative:
//set choose_alt = 0 to generate independent random variables, choose_alt = 1 to generate *epsilon uncorrelated and dependent with *eta1 and *eta2,
//choose_alt = 2 to generate dependent random variables. OUTPUT: three random variables (epsilon, eta1, and eta2).
void NL_Dgp::gen_TriSN (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double rho12, const double rho13, const int choose_alt,
                        unsigned long seed) {
    gsl_rng *r = nullptr;
	const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    double xi_0 = gsl_ran_ugaussian (r);
    Matrix xi(3, 1);
    xi(1) = gsl_ran_ugaussian (r);
    xi(2) = gsl_ran_ugaussian (r);
    xi(3) = gsl_ran_ugaussian (r);
    switch (choose_alt) {
    	case 0: {
    		        double xi_0_star = gsl_ran_ugaussian (r);
    		        epsilon = delta(1) * fabs(xi_0) + sqrt(1 - pow(delta(1),2.)) * xi(1) - sqrt((double) 2/M_PI) * delta(1);
	                eta1    = delta(2) * fabs(xi_0_star) + sqrt(1 - pow(delta(2),2.)) * xi(2) - sqrt((double) 2/M_PI) * delta(2);
	                eta2    = delta(3) * fabs(xi_0_star) + sqrt(1 - pow(delta(3),2.)) * xi(3) - sqrt((double) 2/M_PI) * delta(3);
	                gsl_rng_free (r); //free memory
	            }
    		    break;
    	case 1: {
    		        gsl_matrix *variance = gsl_matrix_alloc (3,3);
    		        gsl_matrix_set (variance, 0, 0, 1.);
    		        double _rho12 = -delta(1)*delta(2)*(1 - 2/M_PI)/sqrt((1 - pow(delta(1),2.))*(1 - pow(delta(2),2.)));
    		        //cout << "_rho12 = " << _rho12 << endl;
    		        gsl_matrix_set (variance, 0, 1, _rho12);
    		        double _rho13 = -delta(1)*delta(3)*(1 - 2/M_PI)/sqrt((1 - pow(delta(1),2.))*(1 - pow(delta(3),2.)));
    		        //cout << "_rho13 = " << _rho13 << endl;
    		        gsl_matrix_set (variance, 0, 2, _rho13);
                    gsl_matrix_set (variance, 1, 0, _rho12);
                    gsl_matrix_set (variance, 1, 1, 1.);
                    gsl_matrix_set (variance, 1, 2, 0.6);
                    gsl_matrix_set (variance, 2, 0, _rho13);
                    gsl_matrix_set (variance, 2, 1, 0.6);
                    gsl_matrix_set (variance, 2, 2, 1.);
                    xi = NL_Dgp::multi_norm (variance, seed); //generate a trivariate normal random variable
                    epsilon = delta(1) * fabs(xi_0) + sqrt(1 - pow(delta(1),2.)) * xi(1) - sqrt((double) 2/M_PI) * delta(1);
	                eta1    = delta(2) * fabs(xi_0) + sqrt(1 - pow(delta(2),2.)) * xi(2) - sqrt((double) 2/M_PI) * delta(2);
	                eta2    = delta(3) * fabs(xi_0) + sqrt(1 - pow(delta(3),2.)) * xi(3) - sqrt((double) 2/M_PI) * delta(3);
	                gsl_matrix_free (variance);
	                gsl_rng_free (r); //free memory
	            }
	            break;
    	case 2: {
                    gsl_matrix *variance = gsl_matrix_alloc (3,3);
                    gsl_matrix_set (variance, 0, 0, 1.);
                    gsl_matrix_set (variance, 0, 1, rho12);
                    gsl_matrix_set (variance, 0, 2, rho13);
                    gsl_matrix_set (variance, 1, 0, rho12);
                    gsl_matrix_set (variance, 1, 1, 1.);
                    gsl_matrix_set (variance, 1, 2, 0.6);
                    gsl_matrix_set (variance, 2, 0, rho13);
                    gsl_matrix_set (variance, 2, 1, 0.6);
                    gsl_matrix_set (variance, 2, 2, 1.);
                    xi = NL_Dgp::multi_norm (variance, seed); //generate a trivariate normal random variable
                    epsilon = delta(1) * fabs(xi_0) + sqrt(1 - pow(delta(1),2.)) * xi(1) - sqrt((double) 2/M_PI) * delta(1);
	                eta1    = delta(2) * fabs(xi_0) + sqrt(1 - pow(delta(2),2.)) * xi(2) - sqrt((double) 2/M_PI) * delta(2);
	                eta2    = delta(3) * fabs(xi_0) + sqrt(1 - pow(delta(3),2.)) * xi(3) - sqrt((double) 2/M_PI) * delta(3);
	                gsl_matrix_free (variance);
	                gsl_rng_free (r); //free memory
	            }
    		    break;
    	default:
    		    cerr << "NL_Dgp::gen_TriSN: This choice is not in the switch list. Make sure that your choice is valid!\n";
                exit(0);
	}
}

//generate two mixtures of N(0,1) random variables. INPUT: a correlation between xi_1_eps and xi_1_eta (rho) and an alternative:
//choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
//dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: two random error terms (epsilon and eta).
void NL_Dgp::gen_MN (double &epsilon, double &eta, const double delta, const double rho, const int choose_alt, unsigned long seed) {
	(void)delta; //unused parameter
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    double xi_0_eps = gsl_ran_ugaussian (r);
    double xi_0_eta = gsl_ran_ugaussian (r);
    double xi_1_eps = gsl_ran_ugaussian (r);
	double xi_1_eta = 0.;
	switch (choose_alt) {
        case 0:
                xi_1_eta = gsl_ran_ugaussian (r); //independent error terms
                break;
        case 1:
        	    gsl_ran_bivariate_gaussian (r, 1., 1., rho, &xi_1_eps, &xi_1_eta); //correlated error terms
        	    break;
        case 2:
        	    if (fabs(xi_1_eps) <= 1.54) //set threshold equal to 1.54 for zero correlation
	                xi_1_eta = xi_1_eps;
	            else
	                xi_1_eta = -xi_1_eps;
	            break;
	    case 3:
	            if (fabs(xi_1_eps) <= 2.) //set threshold equal to 2.0 for non-zero correlation
	                xi_1_eta = xi_1_eps;
	            else
	                xi_1_eta = -xi_1_eps;
	            break;
	    default:
    		    cerr << "NL_Dgp::gen_MN: This choice is not in the switch list. Make sure that your choice is valid!\n";
                exit(0);
	}
	auto z_eps = gsl_ran_bernoulli (r, (double) 1/4);
	if (z_eps == 1)
	    epsilon = xi_1_eps;
	else
	    epsilon = xi_0_eps;
	auto z_eta = gsl_ran_bernoulli (r, (double) 1/3);
	if (z_eta == 1)
	    eta = xi_1_eta;
	else
	    eta = xi_0_eta;
	gsl_rng_free (r); //free memory
}

//generate three mixtures of N(0,1) random variables. INPUT: a correlation between xi_1_eps and xi_1_eta1 (rho) and an alternative:
//choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
//dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: three random error terms (epsilon, eta1, and eta2).
void NL_Dgp::gen_TriMN (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double rho, const double rho1, const int choose_alt,
                        unsigned long seed) {
	(void)delta; //unused parameters
	(void)rho1;
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    double xi_0_eps = gsl_ran_ugaussian (r);
    double xi_0_eta = gsl_ran_ugaussian (r);
    double xi_1_eps = gsl_ran_ugaussian (r);
	double xi_1_eta1 = gsl_ran_ugaussian (r);
	double xi_1_eta2 = gsl_ran_ugaussian (r);
	switch (choose_alt) {
        case 0: //independent error terms
                break;
        case 1:
        	    gsl_ran_bivariate_gaussian (r, 1., 1., rho, &xi_1_eps, &xi_1_eta1); //correlated error terms
        	    break;
        case 2:
        	    if (fabs(xi_1_eps) <= 1.54) {//set threshold equal to 1.54 for zero correlation
	                xi_1_eta1 = xi_1_eps;
	                xi_1_eta2 = xi_1_eps;
	            }
	            else {
	                xi_1_eta1 = -xi_1_eps;
	                xi_1_eta2 = -xi_1_eps;
	            }
	            break;
	    case 3:
	            if (fabs(xi_1_eps) <= 2.0) {//set threshold equal to 2.0 for non-zero correlation
	                xi_1_eta1 = xi_1_eps;
	                xi_1_eta2 = xi_1_eps;
	            }
	            else {
	                xi_1_eta1 = -xi_1_eps;
	                xi_1_eta2 = -xi_1_eps;
	            }
	            break;
	    default:
    		    cerr << "NL_Dgp::gen_TriMN: This choice is not in the switch list. Make sure that your choice is valid!\n";
                exit(0);
	}
	auto z_eps = gsl_ran_bernoulli (r, (double) 1/4);
	if (z_eps == 1)
	    epsilon = xi_1_eps;
	else
	    epsilon = xi_0_eps;
	auto z_eta1 = gsl_ran_bernoulli (r, (double) 1/3);
	if (z_eta1 == 1)
	    eta1 = xi_1_eta1;
	else
	    eta1 = xi_0_eta;
	auto z_eta2 = gsl_ran_bernoulli (r, (double) 2/5);
	if (z_eta2 == 1)
	    eta2 = xi_1_eta2;
	else
	    eta2 = xi_0_eta;
	gsl_rng_free (r); //free memory
}

//generate two centered chi-squared random variables. INPUT: a correlation between xi_1_eps and xi_1_eta (rho) and an alternative:
//choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
//dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: two random error terms (epsilon and eta).
void NL_Dgp::gen_ChiSq (double &epsilon, double &eta, const double delta, const double rho, const int choose_alt, unsigned long seed) {
	(void)delta; //unused parameter
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    double xi_0_eps = gsl_ran_ugaussian (r);
    double xi_0_eta = gsl_ran_ugaussian (r);
    double xi_1_eps = gsl_ran_ugaussian (r);
	double xi_1_eta = 0.;
	switch (choose_alt) {
        case 0:
                xi_1_eta = gsl_ran_ugaussian (r); //independent error terms
                break;
        case 1:
        	    gsl_ran_bivariate_gaussian (r, 1., 1., rho, &xi_1_eps, &xi_1_eta); //correlated error terms
        	    break;
        case 2:
        	    if (fabs(xi_1_eps) <= 1.54) //set threshold equal to 1.54 for zero correlation
	                xi_1_eta = xi_1_eps;
	            else
	                xi_1_eta = -xi_1_eps;
	            break;
	    case 3:
	            if (fabs(xi_1_eps) <= 2.) //set threshold equal to 2.0 for non-zero correlation
	                xi_1_eta = xi_1_eps;
	            else
	                xi_1_eta = -xi_1_eps;
	            break;
	    default:
    		    cerr << "NL_Dgp::gen_ChiSq: This choice is not in the switch list. Make sure that your choice is valid!\n";
                exit(0);
	}
	epsilon = (pow(xi_0_eps, 2.) + pow(xi_1_eps, 2.) - 2)/2;
	eta     = (pow(xi_0_eta, 2.) + pow(xi_1_eta, 2.) - 2)/2;
	gsl_rng_free (r); //free memory
}

//generate three centered chi-squared random variables. INPUT: a correlation between xi_1_eps and xi_1_eta1 (rho) and an alternative:
//choose_alt = 0 to generate independent error terms, choose_alt = 1 to generate correlated error terms, choose_alt = 2 to generate uncorrelated but
//dependent error terms, choose_alt = 3 to generate correlated and dependent error terms. OUTPUT: three random error terms (epsilon, eta1, and eta2).
void NL_Dgp::gen_TriChiSq (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double rho, const double rho1, const int choose_alt,
                           unsigned long seed) {
	(void)delta; //unused parameters
	(void)rho1;
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    double xi_0_eps = gsl_ran_ugaussian (r);
    double xi_0_eta = gsl_ran_ugaussian (r);
    double xi_1_eps = gsl_ran_ugaussian (r);
	double xi_1_eta1 = gsl_ran_ugaussian (r);
	double xi_1_eta2 = gsl_ran_ugaussian (r);
	switch (choose_alt) {
        case 0: //independent error terms
                break;
        case 1:
        	    gsl_ran_bivariate_gaussian (r, 1., 1., rho, &xi_1_eps, &xi_1_eta1); //correlated error terms
        	    break;
        case 2:
        	    if (fabs(xi_1_eps) <= 1.54) {//set threshold equal to 1.54 for zero correlation
	                xi_1_eta1 = xi_1_eps;
	                xi_1_eta2 = xi_1_eps;
	            }
	            else {
	                xi_1_eta1 = -xi_1_eps;
	                xi_1_eta2 = -xi_1_eps;
	            }
	            break;
	    case 3:
	            if (fabs(xi_1_eps) <= 2.0) {//set threshold equal to 2.0 for non-zero correlation
	                xi_1_eta1 = xi_1_eps;
	                xi_1_eta2 = xi_1_eps;
	            }
	            else {
	                xi_1_eta1 = -xi_1_eps;
	                xi_1_eta2 = -xi_1_eps;
	            }
	            break;
	    default:
    		    cerr << "NL_Dgp::gen_TriMN: This choice is not in the switch list. Make sure that your choice is valid!\n";
                exit(0);
	}
	epsilon = (pow(xi_0_eps, 2.) + pow(xi_1_eps, 2.) - 2)/2;
	eta1     = (pow(xi_0_eta, 2.) + pow(xi_1_eta1, 2.) - 2)/2;
	eta2     = (pow(xi_0_eta, 2.) + pow(xi_1_eta2, 2.) - 2)/2;
	gsl_rng_free (r); //free memory
}

//generate two centered Beta random variables. INPUT: a shape parameter (alpha5) -- alpha5 = 0.0001, 0.1, 1., 2. OUTPUT: two random error terms (epsilon and eta).
void NL_Dgp::gen_Beta (double &epsilon, double &eta, const double alpha5, const double rho, const int choose_alt, unsigned long seed) {
	(void)rho; //unused parameters
	(void)choose_alt;
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    Matrix U(5,1);
    U(1) = gsl_ran_gamma (r, 5., 1.);
    U(2) = gsl_ran_gamma (r, 5., 1.);
    U(3) = gsl_ran_gamma (r, 0.5, 1.);
    U(4) = gsl_ran_gamma (r, 0.5, 1.);
    U(5) = gsl_ran_gamma (r, alpha5, 1.);
    epsilon = ((double) (U(1) + U(3))/(U(1) + U(3) + U(4) + U(5)) - (5 + 0.5)/(5 + 0.5 + 0.5 + alpha5));
    eta     = ((double) (U(2) + U(4))/(U(2) + U(3) + U(4) + U(5)) - (5 + 0.5)/(5 + 0.5 + 0.5 + alpha5));
    gsl_rng_free (r);//free memory
}

//generate three centered Beta random variables. INPUT: a shape parameter (alpha5) -- alpha5 = 0.0001, 0.1, 1., 2.
//OUTPUT: three random error terms (epsilon, eta1, and eta2).
void NL_Dgp::gen_TriBeta (double &epsilon, double &eta1, double &eta2, const Matrix &delta, const double alpha5, const double rho1, const int choose_alt,
                          unsigned long seed) {
	(void)delta; //unused parameters
	(void)rho1;
	(void)choose_alt;
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    Matrix U(5,1);
    U(1) = gsl_ran_gamma (r, 5., 1.);
    U(2) = gsl_ran_gamma (r, 5., 1.);
    U(3) = gsl_ran_gamma (r, 0.5, 1.);
    U(4) = gsl_ran_gamma (r, 0.5, 1.);
    U(5) = gsl_ran_gamma (r, alpha5, 1.);
    epsilon = ((double) (U(1) + U(3))/(U(1) + U(3) + U(4) + U(5)) - (5 + 0.5)/(5 + 0.5 + 0.5 + alpha5));
    eta1     = ((double) (U(2) + U(4))/(U(2) + U(3) + U(4) + U(5)) - (5 + 0.5)/(5 + 0.5 + 0.5 + alpha5));
    eta2 = 0.5*eta1 + gsl_ran_ugaussian (r);
    gsl_rng_free (r);//free memory
}

//generate two vectors of random errors (epsilon and eta). INPUT: constants (delta and rho in (-1,1)), an alternative (choose_alt = 0, 1, 2), and a
//random generator template (gen_RAN: gen_SN, gen_MN, gen_ChiSq or gen_Beta). OUTPUT: 2 vectors
template <void gen_RAN (double &, double &, const double, const double, const int, unsigned long)>
void NL_Dgp::gen_RANV (Matrix &epsilon, Matrix &eta, const double delta, const double rho, const int choose_alt, unsigned long seed) {
	gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    unsigned long rseed = 1;
    for (auto i = 1; i <= epsilon.nRow(); i++) {
    	rseed = gsl_rng_get (r);
    	gen_RAN (epsilon(i), eta(i), delta, rho, choose_alt, rseed); //use a random seed
	}
    gsl_rng_free (r); //free memory
}

//generate three vectors of random errors (epsilon, eta1 and eta2). INPUT: a vector (delta in (-1,1)^3) and constants (rho12 and rho13 in (-1,1)),
//an alternative (choose_alt = 0, 1, 2), and a random generator template (gen_RAN: gen_TriSN, gen_TriMN, gen_TriChiSq or gen_TriBeta). OUTPUT: 3 vectors
template <void gen_RAN (double &, double &, double &, const Matrix &, const double, const double, const int, unsigned long)>
void NL_Dgp::gen_RANV (Matrix &epsilon, Matrix &eta1, Matrix &eta2, const Matrix &delta, const double rho12, const double rho13, const int choose_alt,
                       unsigned long seed) {
    gsl_rng *r = nullptr;
    const gsl_rng_type *gen; //random number generator
    gsl_rng_env_setup();
    gen = gsl_rng_taus;
    r = gsl_rng_alloc(gen);
    gsl_rng_set(r, seed);
    unsigned long rseed = 1;
    for (auto i = 1; i <= epsilon.nRow(); i++) {
    	rseed = gsl_rng_get (r);
    	gen_RAN (epsilon(i), eta1(i), eta2(i), delta, rho12, rho13, choose_alt, rseed);
    	//cout << epsilon (i) << " , " << eta1(i) << " , " << eta2(i) << endl;
	}
    gsl_rng_free (r); //free memory

}

//run multivariate OLS. INPUT: a Tx1 vector of data on the dependent (Y) and a TxN matrix of data on the independent (X).
//OUTPUT: a Tx1 vector of residuals (resid) and a Nx1 vector of the OLS estimates (slope)
void NL_Dgp::gen_Resid (Matrix &resid, Matrix &slope, const Matrix X, const Matrix Y) {
	//auto T = X.nRow();
	auto N = X.nCol();
	Matrix denom(N,N), num(N,1), denom_Inv(N,N);
	denom = Tr(X)*X;
	denom_Inv = inv(denom);
	num = Tr(X)*Y;
	slope = denom_Inv*num;
	resid = Y - X*slope;
}

//estimate the bilinear regression model by the OLS. INPUT: a Tx1 vector of data on X.
//OUTPUT: a (T-2)x1 vector of residuals (resid) and a 6x1 vector of the OLS estimates (slope)
void NL_Dgp::est_BL (Matrix &resid, Matrix &slope, const Matrix X) {
	auto T = X.nRow(), row = T-2, col = 6;
	Matrix Z(row,col), Y(row,1);
	for (auto t = 1; t <= row; ++t) {
		Z(t,1) = 1.;
		Z(t,2) = X(t+1);
		Z(t,3) = pow(Z(t,2), 2.);
		Z(t,4) = X(t);
		Z(t,5) = pow(Z(t,4), 2.);
		Z(t,6) = Z(t,2) * Z(t,4);
		Y(t) = X(t+2);
	}
	NL_Dgp::gen_Resid (resid, slope, Z, Y);
}

//estimate the TAR regression model by the OLS. INPUT: a Tx1 vector of data on X.
//OUTPUT: a (T-2)x1 vector of residuals (resid) and a 5x1 vector of the OLS estimates (slope)
void NL_Dgp::est_TAR (Matrix &resid, Matrix &slope, const Matrix X) {
	auto T = X.nRow(), row = T-2, col = 5;
	Matrix Z(row,col), Y(row,1);
	auto pos = [](double x) { //get the positive part of a real number (x)
    	if (x > 0.)
    		return x;
    	else
    		return 0.;
	};
	for (auto t = 1; t <= row; ++t) {
		Z(t,1) = 1.;
		Z(t,2) = X(t+1);
		Z(t,3) = pos(Z(t,2));
		Z(t,4) = X(t);
		Z(t,5) = pos(Z(t,4));
		Y(t) = X(t+2);
	}
	NL_Dgp::gen_Resid (resid, slope, Z, Y);
}






















#endif