testreg.cpp

#include "virt.h"

//required global declaration for the
//  matrix stack object

unsigned int _stklen = STACKLENGTH;
MStack *Dispatch = new MStack;

VMatrix& ExactHilbInv( int n )
{
   // construct exact hilbert matrix inverse
   Dispatch->Inclevel();
   VMatrix Hi("Hilbert Inverse",n,n);

   double diag, dn = n;
   for( int i=1; i<=n; i++){
	  double im1 = (double) (i-1.0);
	  diag = ( i == 1 ) ? dn : ((dn-im1)*diag*(dn+im1))/(im1*im1);

	  long double rest = diag*diag;
	  Hi.M(i,i) = rest/(2.0*im1+1.0);
	  for( int j=i+1; j<=n; j++ ){
		 double jm1 = (double) (j-1);
		 rest = -((dn-jm1)*rest*(dn+jm1))/(jm1*jm1);
		 Hi.M(i,j) = rest/(im1+jm1+1.0);
		 Hi.M(j,i) = Hi.m(i,j);
	  }
   }
   Dispatch->Push( Hi );
   return Dispatch->DecReturn();
}


void TestHilb( int hilb_ind )
// test svd on hilbert matrix
{
  Dispatch->Inclevel();
  VMatrix hilb = Kron( Fill(1,hilb_ind,1), Index(1, hilb_ind));
  // h(i,j) = 1.0/(i+j-1)
  hilb = 1.0 / ( (hilb+Tran(hilb))-1.0);
  hilb.Nameit("Hilbert Matrix");
  hilb.DisplayMat();
  // use formula. Looses accuracy of 0.001 at n>=11
  (ExactHilbInv( hilb_ind )*hilb ).DisplayMat();

  // use sweep
  VMatrix Gauss = Inv(hilb)*hilb;
  Gauss.Nameit("Gaussian elimination");
  Gauss.DisplayMat();

  // use singular valued decomposition
  VMatrix S, V, D;
  Svd( hilb, S, V, D);
  VMatrix t = Fill(V.r, V.r, 0.0);
  for (int i = 1; i <= V.r; i++) {
	double tt = V.m(i, 1);
	t.M(i, i) = 1.0 / tt;
  }

  VMatrix Ggauss = (D*t*Tran(S))*hilb;
  Ggauss.Nameit("Generalized inversion");
  Ggauss.DisplayMat();

  Dispatch->Declevel();
}

///////////////////// now test regressions on sample data

VMatrix &getx(int N)
// create an x matrix
{
	Dispatch->Inclevel();

	VMatrix x = Index( 1, N ) - N/2;
	VMatrix c1 = Fill(N, 1, 1.0);
	VMatrix x2 = x % x;
	x = Ch(c1, Ch(x, x2));

	// push x onto the stack
	Dispatch->Push(x);
	// decrement the subroutine nesting level
	// and return the stack top
	return Dispatch->DecReturn();
}

VMatrix &gety(VMatrix &x, VMatrix &beta)
// create a y matrix
{
	Dispatch->Inclevel();

	VMatrix y = x*beta;
	srand(123);
	for (int i = 1; i <= y.r; i++) {
		// use sum of 3 uniforms for an approximate
		// normal random variable
		int u = random(100) + random(100) + random(100) + 3;
		y.M(i, 1) = y.m(i, 1) + ((double) (u - 150)) / 300.0;
	}
	Dispatch->Push(y);
	return Dispatch->DecReturn();
}

VMatrix &regression(VMatrix& x, VMatrix& y)
// do a multiple linear regression
{
	Dispatch->Inclevel();

	int N = x.r, p = x.c;
	// solve for regression parameters using sweep
	VMatrix yx = Ch(y, x);
	VMatrix reg = Sweep(2, p + 1, Tran(yx) *yx);
	// calculate mean square error
	reg.M(1, 1) = reg.m(1, 1) / ((double) (N - p));
	reg.DisplayMat();


	// solve regression using normal equations
	VMatrix betahat = Inv(Tran(x) *x) *Tran(x) *y;

	Dispatch->Push(betahat);
	return Dispatch->DecReturn();
}
VMatrix &QRregression(VMatrix &x, VMatrix &y)
{
	// use QR decomposition to solve regression

	Dispatch->Inclevel();
	VMatrix Q, R;

	QR(x, Q, R);
	VMatrix betahat = Inv(Tran(R) *R) *Tran(R) *Tran(Q) *y;

	Dispatch->Push(betahat);
	return Dispatch->DecReturn();
}

VMatrix &GinvRegression(VMatrix &x, VMatrix &y)
{
	// use Ginv to solve normal equations
	Dispatch->Inclevel();

	VMatrix betahat = Ginv(Tran(x) *x) *Tran(x) *y;

	Dispatch->Push(betahat);
	return Dispatch->DecReturn();
}


VMatrix &SVDregression(VMatrix &x, VMatrix &y)
{
	// use SVD to solve regression x = S Diag(V) Tran(D)
	Dispatch->Inclevel();
	VMatrix S, V, D;

	Svd(x, S, V, D);
	VMatrix t = Fill(V.r, V.r, 0.0);
	for (int i = 1; i <= V.r; i++) {
		double tt = V.m(i, 1);
		t.M(i, i) = (fabs(tt) <= 0.0) ? 0.0 : 1.0 / tt;
	}
	VMatrix betahat = D*t*Tran(S)*y;

	Dispatch->Push(betahat);
	return Dispatch->DecReturn();
}
VMatrix &GetSerialCovar(VMatrix &R)
{
   // Parameters to CORR in Timeslab
   // correlogram = CORR(x=R,n=R.r,M=R.r-1,Q=2*R.r,
   //                    iopt=1,r0=r0, per = spectdens)
   // covar = correlogram*r0
   Dispatch->Inclevel();
   VMatrix centered, z, covar, spectdens;
   double n = (double) R.r;
   // center a column vector
   centered = R - Sum(R).m(1, 1) / n;
   // zero pad to length 2n: 2n periodic for full
   // sample spectral density
   centered = Cv(centered, Fill(R.r, R.c, 0));
   // take fft
   z = Fft(centered);
   // take convolution : gives sample spectral density
   spectdens = Sum(z % z / n, COLUMNS);
   // inverse fft for serial correlation (autocorrelation)
   covar = Fft(spectdens, FFALSE);
   // throw away last half.
   covar = Submat(covar, R.r, 2);
   Dispatch->Push(covar);
   return Dispatch->DecReturn();
}

main()
{
	Dispatch->Inclevel();
	VMatrix x, beta("beta", 3, 1), y, betahat, resids, serial;

	TestHilb( 11 );

	Setwid(15);
	Setdec(10);


	beta.M(1, 1) = 1;
	beta.M(2, 1) = 0.5;
	beta.M(3, 1) = 0.25;

	x = getx(55);
	y = gety(x, beta);

	betahat = regression(x, y);
	betahat.Nameit("Text book betahat");
	(Tran(beta)).DisplayMat();
	(Tran(betahat)).DisplayMat();

	betahat = QRregression(x, y);
	betahat.Nameit("QR betahat");
	(Tran(beta)).DisplayMat();
	(Tran(betahat)).DisplayMat();

	betahat = GinvRegression(x, y);
	betahat.Nameit("Ginv betahat");
	(Tran(beta)).DisplayMat();
	(Tran(betahat)).DisplayMat();

	betahat = SVDregression(x, y);
	betahat.Nameit("SVD regression");
	(Tran(beta)).DisplayMat();
	(Tran(betahat)).DisplayMat();

	resids = y - x*betahat;
	serial = GetSerialCovar(resids);
	(Tran(Submat(serial, 5, 1))).DisplayMat();

	Setwid(6);
	Setdec(3);

	vclose();
}