\chapter{testgraf.cpp}
This is the test suite for the graphics routines. It plots the
residuals, and correlogram and periodogram of the residuals. See
Newton \cite{NJ:tl} or Press et.al. \cite{Pr:nrc} 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. 

\begin{verbatim}


#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%4))/3.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 )
  {
    Dispatch->Inclevel();
    VMatrix centered, z, covar;
    double n = (double) R.r;
    centered = R - Sum( R/n ).m(1,1);
    z = Fft( centered );
    spectdens = Sum( z%z/n, COLUMNS);
    covar = Fft( spectdens, FFALSE );
    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/2, 1, 2, 1);
     VMatrix graf = Ch( Index( 1, Correlogram.r ), Correlogram );

     GMatrix Agraph( graf );
     double dn2 = 2.0/sqrt((double)n);
     *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 )
  {
     Dispatch->Inclevel();
     int n = spectdens.r;
     double dn = (double)n;
     double sigma = (Sum(spectdens).m(1,1) - spectdens.m(1,1));
     VMatrix Periodogram = spectdens/sigma;
     Periodogram = Mlog(Submat( Periodogram, 1+n/2, 1, 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();
}

\end{verbatim}