cppcode.cpp

#include <random>
#include <complex>
#include <iostream>
#include <thrust/host_vector.h>
#include "armadillo"
#include <lapacke.h>

using namespace arma;
using namespace std;

void mynormalcpp(double* output, double mmu, double ssigma, int n, unsigned seed) {
	
	// Specify engine and distribution
	std::default_random_engine generator(seed);
	std::normal_distribution<double> distribution(mmu,ssigma);

	// Actually generate random numbers
	for (int i=0; i<n; i++) {	
		output[i] = distribution(generator);
	};
	
};

void myuniformcpp(double* output, int n, unsigned seed)
{
	std::default_random_engine generator(seed);
    std::uniform_real_distribution<double> d(0,1);
    
	for (int i=0; i<n; i++) {
		output[i] = d(generator);
	};
};

void myexponentialcpp(double* output, int n, unsigned seed)
{
	std::default_random_engine generator(seed);
    std::exponential_distribution<double> d(1.0);
    
	for (int i=0; i<n; i++) {
		output[i] = d(generator);
	};
};

int chebyroots_cpp(const int p, double* roots) {
    for (int i=0; i<p; i++) {
        double stuff = p - 0.5 - 1*i;
        roots[i] = cos(M_PI*(stuff)/(p));
    };

    // Account for the fact that cos(pi/2) is not exactly zeros
    if (p%2) {
        roots[(p-1)/2] = 0;
    };
    return 0;
};

void fromchebydomain(double lb, double ub, int p, double* ptr) {
	for (int i=0; i<p; i++) {
		ptr[i] = lb + (ptr[i]+1)/(2)*(ub-lb);	
	};
};

double mynormalcdf(double mean, double variance, double x) {
	return 0.5*( 1+erf( (x-mean)/sqrt(2*variance) ) );
};

// Create a equal-distance grid
void linspace(double min, double max, int N, double* grid) {
	double step  = (max-min)/(N-1);
	for (int i = 0; i < N; i++) {
		grid[i] = min + step*i;
	};
};

// Create a grid based on chebyshev roots
void chebyspace(double min, double max, int N, double* grid) {
	chebyroots_cpp(N, grid);
	fromchebydomain(min, max, N, grid);
};

// Define a type of function pointer
typedef void (*gridgen_fptr)(double, double, int, double*);

// Declare ind2sub here. Code is generated by cuda_helpers.h and compiled by nvcc though. The linker will find the right machine code for me. Now this code is dependent on cuda_helpers.h
void ind2sub(int, int*, int, int*);

// This function follows notation from Tauchen 1986 the vector case using Armadillo
void tauchen_vec(int M, int N, int m, double* A_ptr, double* Ssigma_e_ptr, double* Z_ptr, double* P_ptr, gridgen_fptr gridgen) {
// Purpose:		Discretize y_t = A*y_t-1 + epsilon_t, where var(epsilon_t)=Ssigma_e
// 				y_t is an M by 1 vector
//
// Input:		M = # of shocks
// 				N = # of grid points for each shock
// 				m = # of s.d. away from mean should we use as bounds 
// 				A = Autocorrelation matrix, stored in array 
// 				Ssigma_e = variance-convariance matrix, stored in array
//
// Output:		Z = Vectorized N by M matrix where each column stores grids
// 				P = Vectorized transition matrix (N^M by N^M)
	
	// 1. Read in matrices
	mat A(M,M);
	mat Ssigma_e(M,M);
	for (int i=0; i < M*M; i++) {
		int i_c = i/M;
		int i_r = i - i_c*M;
		A(i_r,i_c) = A_ptr[i];
		Ssigma_e(i_r,i_c) = Ssigma_e_ptr[i];
	};

	// 2. Find Variance Matrix for y
	vec Ssigma_e_vec = vectorise(Ssigma_e);
	vec Ssigma_y_vec = inv( eye(M*M,M*M) - kron(A,A) ) * Ssigma_e_vec;
	mat Ssigma_y = reshape(Ssigma_y_vec,M,M);
	A.print("Autocorr Matrix:");
	Ssigma_e.print("Var-Cove Matrix of Innovations:");
	Ssigma_y.print("Var-Cov Matrix of Observables:");

	// 3. Create Grids. Assuming same # of grid points for each shock
	mat grids(M,N); // collects grid points here
	for (int i_shock=0; i_shock<M; i_shock++) {
		double minshock = -m*sqrt(Ssigma_y(i_shock,i_shock));
		double maxshock = +m*sqrt(Ssigma_y(i_shock,i_shock));
		vec temp_grid(N);
		gridgen(minshock,maxshock,N,temp_grid.memptr());
		grids.row(i_shock) = trans(temp_grid);
	};
	grids.print("Grids:");

	// 4. Compute Transition Matrix
	mat P(pow(N,M),pow(N,M));
	cube h(M,pow(N,M),N);
	Col<int> sizes(M); sizes.fill(N);
	Col<int> subs(M);
	vec lag_y(M);
	for (int j=0; j<pow(N,M); j++) {
		// First find the conditional mean
		ind2sub(M, sizes.memptr(), j, subs.memptr());
		for (int i_shock=0; i_shock<M; i_shock++) lag_y(i_shock) = grids(i_shock,subs(i_shock));
		vec mmu = A*lag_y;

		// Fill the h cube
		double w_right, w_left;
		for (int i=0; i<M; i++) {
			for (int l=0; l<N; l++) {
				if (l==0) {
					w_right = grids(i,l+1) - grids(i,l);
					h(i,j,l) = mynormalcdf(0,Ssigma_e(i,i),grids(i,l)-mmu(i)+w_right/2);
				}
			   	else if (l==N-1) {
					w_left = grids(i,l) - grids(i,l-1);
					h(i,j,l) = 1 - mynormalcdf(0,Ssigma_e(i,i),grids(i,l)-mmu(i)-w_left/2);
				}
			   	else {
					w_left = grids(i,l) - grids(i,l-1);
					w_right = grids(i,l+1) - grids(i,l);
					h(i,j,l) = mynormalcdf(0,Ssigma_e(i,i),grids(i,l)-mmu(i)+w_right/2)
							 - mynormalcdf(0,Ssigma_e(i,i),grids(i,l)-mmu(i)-w_left/2);
				};
			};
		};
	};
	for (int j=0; j<pow(N,M); j++) {
		for (int k=0; k<pow(N,M); k++) {
			ind2sub(M, sizes.memptr(), k, subs.memptr());
			double cumprod = 1;
			for (int i=0; i<M; i++) {
				cumprod *= h(i,j,subs(i));
			};
			P(j,k) = cumprod;
		};
	};
	vec rowsum = sum(P,1);
	
	// 5. Check Accuracy
	arma_rng::set_seed(51709394);
	int T = 1e3;
	mat eps(M,T);
	for (int i_shock=0; i_shock<M; i_shock++) { 
		eps(i_shock,span::all) = trans(sqrt(Ssigma_e(i_shock,i_shock))*randn<vec>(T));
	};
	mat sim_y = zeros<mat>(M,T);
	for (int t=0; t<T-1; t++) {
		sim_y(span::all,t+1) = A*sim_y(span::all,t) + eps(span::all,t+1);
	};
	mat Z = sim_y(span::all,span(0,T-2));
	mat Y = sim_y(span::all,span(1,T-1));
	mat B(M,M);
	B = Y*trans(Z)*inv(Z*trans(Z)); //Multivariate Regression, formulas from Wikipedia
	B.print("Discrete Autocorr Matrix:");
	mat C = cov(trans(sim_y));
	C.print("Discrete Var-Cov Matrix of Observables: ");

	// 6. Output Results
	mat Z_out = trans(grids);
	double *Z_out_ptr = Z_out.memptr();
	for (int index=0; index<N*M; index++) {
		Z_ptr[index] = Z_out_ptr[index];
	};
	double* P_out_ptr = P.memptr();
	for (int index=0; index<pow(N,M)*pow(N,M); index++) {
		P_ptr[index] = P_out_ptr[index]; 
	};
};

void findprojector(double* X_ptr, int nrows, int ncols, double * P_ptr) {
	mat X(X_ptr, nrows, ncols, false, false);
	mat XprimeX = diagmat(trans(X)*X);
	mat P = inv(XprimeX)*trans(X);
	
	for (int i=0; i<nrows*ncols; i++) {
		P_ptr[i] = P.memptr()[i];
	};

};

lapack_logical qzdivide(const complex<double>* alpha, const complex<double>* beta) {
	// As in Sims and Klein, our Generalised eigenvalues are defined as b/a 
	complex<double> a = * alpha;
	complex<double> b = * beta;
	complex<double> zero(0,0);

	// Determined whether eigenvalues are within unit circle, assuming not all a,b,c are zeros.
	// We want the >1 on upper left, because we will invert it later.
	if (a==zero) {
		return false;
	} else if (abs(b/a)<1) {
		return true;
	} else {
		return false;
	};
};

int qzdecomp (cx_mat &A, cx_mat &B, cx_mat &Q, cx_mat &Z) {
// Purpose:
// 	Compute QZ decomposition for a pair nonsymmetric matrice (A,B).
// 	A = Q'*S*Z, B = Q'*T*Z
// 	where S, T is stored in A, B when done.
// 	The decomposition is ordered such that explosive eigenvalues are placed in lower right of inv(S)*T;
// 	Returns the number of unstable roots

	// Preparations
	int n_selected;
	cx_vec alpha(A.n_rows);
	cx_vec beta(B.n_rows);
	lapack_int info;	// stores the exit information 
	LAPACK_Z_SELECT2 selctg = qzdivide; 
	A.print("Inputted A matrix:");
	B.print("Inputted B matrix:");

	// Call LAPACK
	info = LAPACKE_zgges(
			LAPACK_COL_MAJOR,	// indicate we use column major
			'V',				// ... and compute left schur vector 
			'V',				// ... and compute right schur vector
			'S',				// ... and sort eigenvalues
			selctg,				// with this func to divide eigenvalues
			A.n_rows,			// # of rows in A
			A.memptr(),			// array that stores A
			A.n_rows,			// leading dimension of A (rows)
			B.memptr(),			// array that stores B
			B.n_rows,			// rows of B
			&n_selected,		// output the # of eigenvalues 
			alpha.memptr(),	// stores aux vectors
			beta.memptr(),		// stores aux vector
			Q.memptr(),			// stores Q
			A.n_rows,
			Z.memptr(),
			B.n_rows
			);
	
	// Transpose Q consistent with Klein's notation.
	Q = trans(Q);

	// Print and See
	A.print("S = ");
	B.print("T = ");
	Q.print("Q = ");
	Z.print("Z = ");
	printf("The %i ordered eigenvalues have the following norms: \n", A.n_rows);
	complex<double> zero(0,0);
	for (int i = 0; i<A.n_rows; i++) {
		if (beta(i) != zero) {
			printf("[#%i] eigenvalue has norm: %f\n", i+1, abs(alpha(i)/beta(i)) );
		} else {
			printf("#[%i] eigenvalue has norm: Inf\n", i+1  );
		};

	};

	// Check accuray
	cx_mat Atilde = trans(Q)*A*trans(Z);
	Atilde.print("Atilde = ");
	cx_mat Btilde = trans(Q)*B*trans(Z);
	Btilde.print("Btilde = ");

	
	// Compute the number of unstable roots (in lower right of inv(S)*T)
	return A.n_rows - n_selected;
};

void linearQZ(double* A_ptr, double* B_ptr, double* C_ptr, double* rrho_ptr, int n, int n_jump, int n_shock, double* Pphi_ptr) {
// Purpose: (Following Juan's notation)
// 	Solves any linear rational expectation model in the form of:
// 	A* E_t y_t+1 = B* y_t + C * z_t and
// 	z_t = rrho* z_t-1 + G * epsilon_t
// 	y_t is the n-by-1 vector of endogenuous variables which contains n_predeter
// 	predetermined (jump) variables. A, B are n-by-n. C is n-by-n_shock, rrho is n_shock-by-n_shock.
// 	Pphi is n-by-(n-n_jump+n_shock)

	// Wrap armadillo mat around raw arrays and convert to complex matrices.
	mat A_real(A_ptr,n,n) ;
	mat B_real(B_ptr,n,n) ;
	mat C_real(C_ptr,n,n_shock) ;
	mat rrho_real(rrho_ptr,n_shock,n_shock) ;
	cx_mat A = conv_to<cx_mat>::from(A_real); 
	cx_mat B = conv_to<cx_mat>::from(B_real); 
	cx_mat C = conv_to<cx_mat>::from(C_real); 
	cx_mat rrho = conv_to<cx_mat>::from(rrho_real); 
	
	// Call QZ decomposition
	cx_mat S = A;
	cx_mat T = B;
	cx_mat Q(n,n);
	cx_mat Z(n,n);
	int n_unstable = qzdecomp(S,T,Q,Z);

	// Check counting rule
	if (n_unstable > n_jump) {
		printf("There is %i unstables roots but %i jump variables => no solution.\n", n_unstable, n_jump);
	} else if (n_unstable < n_jump) {
		printf("There is %i unstables roots but %i jump variables => infinitely many solutions.\n", n_unstable, n_jump);
	} else {
		printf("There is %i unstables roots but %i jump variables => unique solution.\n", n_unstable, n_jump);
	};

	// Find solution.
	cx_mat S11 = S(span(0,n-n_jump-1),span(0,n-n_jump-1));
	cx_mat S12 = S(span(0,n-n_jump-1),span(n-n_jump,n-1));
	cx_mat S21 = S(span(n-n_jump,n-1),span(0,n-n_jump-1));
	cx_mat S22 = S(span(n-n_jump,n-1),span(n-n_jump,n-1));
	
	cx_mat T11 = T(span(0,n-n_jump-1),span(0,n-n_jump-1));
	cx_mat T12 = T(span(0,n-n_jump-1),span(n-n_jump,n-1));
	cx_mat T21 = T(span(n-n_jump,n-1),span(0,n-n_jump-1));
	cx_mat T22 = T(span(n-n_jump,n-1),span(n-n_jump,n-1));

	cx_mat Z11 = Z(span(0,n-n_jump-1),span(0,n-n_jump-1));
	cx_mat Z12 = Z(span(0,n-n_jump-1),span(n-n_jump,n-1));
	cx_mat Z21 = Z(span(n-n_jump,n-1),span(0,n-n_jump-1));
	cx_mat Z22 = Z(span(n-n_jump,n-1),span(n-n_jump,n-1));

	cx_mat Q1 = Q(span(0,n-n_jump-1),span::all);
	cx_mat Q2 = Q(span(n-n_jump,n-1),span::all);

	// From now on just copying formulas in Klein 2000
	cx_vec M_vec = inv( kron(trans(rrho),S22)-kron(eye(n_shock,n_shock),T22) ) *vectorise(Q2*C);
	cx_mat M = reshape(M_vec,n_jump,n_shock);
	cx_mat N = (Z22-Z21*inv(Z11)*Z12)*M;
	cx_mat L = -Z11*inv(S11)*T11*inv(Z11)*Z12*M + Z11*inv(S11)*(T12*M-S12*M*rrho+Q1*C) + Z12*M*rrho;
	cx_mat Pphi_c_k = Z21*inv(Z11);
	cx_mat Pphi_c_z = N; 
	cx_mat Pphi_k_k = Z11*inv(S11)*T11*inv(Z11); 
	cx_mat Pphi_k_z = L; 
	cx_mat Pphi_cx = join_cols( join_rows(Pphi_k_k,Pphi_k_z), join_rows(Pphi_c_k,Pphi_c_z)  );
	mat Pphi = conv_to<mat>::from(Pphi_cx);
	Pphi.print("The solution matrix is: ");

	// Output Pphi, the solution.
	for (int i = 0; i < n*(n-n_jump+n_shock); i++) {
		Pphi_ptr[i] = Pphi.memptr()[i];
	};
};