testgraf.cpp

This is the test suite for the graphics routines. It plots the residuals, and correlogram and periodogram of the residuals. See Newton [20] or Press et.al. [21] for an explanation of these calculations. They identify serial correlation in the residuals of regression that remove a linear trend from the data. The serial corelation was intentionally placed on the residuals.

#include "virt.h"

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

// unsigned int _stklen = STACKLENGTH;

MStack *Dispatch = new MStack;


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

      c1 = Fill(N,1,1.0);
      x = Index( 1, N ) - ((double)N)*0.5;
      x = Ch( c1,x );

      // 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;

      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) + 5.0*sin( 3.14*((double) (i%8))/7.0 )
                 + ((double) (u-150))/300.0;
      }
      Dispatch->Push(y);
      return Dispatch->DecReturn();
  }

VMatrix &regression( VMatrix& x, VMatrix& y )
// do a multiple linear regression
      {
      Dispatch->Inclevel();
      VMatrix yx, reg, betahat;
      int N=x.r, p=x.c;

      // solve for regression parameters using sweep
      yx = Ch(y,x);
      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
      betahat = Inv(Tran(x)*x)*Tran(x)*y;

      Dispatch->Push( betahat );
      return Dispatch->DecReturn();
}
void plotResiduals( VMatrix &resids )
  {
    Dispatch->Inclevel();
    VMatrix grf = Ch( Index( 1, resids.r ), resids );
    GMatrix Agraph( grf, -'%' );
    *Agraph.PathToDriver = "C:\\tc\\bgi";
    *Agraph.title = "Residuals for data";
    *Agraph.title2= "with serial correlations with frequency 0.25";
    *Agraph.yname = "Residuals";
    *Agraph.xname = "Index";
    Agraph.Href( 0.0 );
    Agraph.Show();
    Dispatch->Cleanstack();
    Dispatch->Declevel();
  }

VMatrix &GetSerialCovar(VMatrix &R, VMatrix &spectdens)
{
   // 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;
   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();
}

void plotCorrelogram(VMatrix &serial)
{
   Dispatch->Inclevel();
   
   int n = serial.r;
   double sigma = serial.m(1, 1);
   VMatrix Correlogram = serial / sigma;
   Correlogram = Submat(Correlogram, n, 1, 2, 1);
   VMatrix graf = Ch(Index(1, Correlogram.r), Correlogram);
   
   GMatrix Agraph(graf);
   *Agraph.PathToDriver = "C:\\tc\\bgi";
   *Agraph.title = "Serial Correlations";
   *Agraph.title2 = "for sample residuals";
   *Agraph.yname = "Serial correlations";
   *Agraph.xname = "Lags";
   Agraph.Href(0.0);
   Agraph.Show();
   
   Dispatch->Cleanstack();
   Dispatch->Declevel();
}

void plotPeriodogram(VMatrix &spectdens)
{
   // calculate a standardized periodogram on the log scale
   Dispatch->Inclevel();
   int n = spectdens.r;
   // this works because data is already centered, which
   // forces spectdens.m(1,1) = 0.0;
   double sigma = Sum(spectdens).m(1, 1) / ((double) n);
   VMatrix Periodogram = spectdens / sigma;
   Periodogram = Mlog(Submat(Periodogram, n / 2 + 1, 1, 2, 1));
   
   // frequencies
   double dn =(double) (n / 2 + 1);
   VMatrix graf = Ch(Index(1, Periodogram.r) / dn, Periodogram);
   
   GMatrix Agraph(graf);
   *Agraph.PathToDriver = "C:\\tc\\bgi";
   *Agraph.title = "Periodogram";
   *Agraph.title2 = "for sample residuals";
   *Agraph.yname = "Log periodogram";
   *Agraph.xname = "Frequencies";
   Agraph.Vref(0.25);
   Agraph.Show();
   
   Dispatch->Cleanstack();
   Dispatch->Declevel();
}


main()
{
      Dispatch->Inclevel();
      VMatrix x, beta("beta",2,1), y, betahat, resids, serial;
      Setwid(15);
      Setdec(10);

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

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

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

      resids = y - x*betahat;
      plotResiduals( resids );

      VMatrix spectdens;
      serial = GetSerialCovar( resids, spectdens );

      plotCorrelogram( serial );
      plotPeriodogram( spectdens);

   vclose();
}