/****************   PVALPLOT       *********************
This procedure I've written that might be useful to others
doing Monte Carlo experiments. It accepts a matrix of test
statistics and information about their theoretical distribution under
null.  It then produces graphs of the extent of size distortion in the
simulated statistics, and returns some descriptive statistics.

The graphs are the p-value and the p-value discrepancy plots suggested
by Davidson and MacKinnon (1997).
    available at: http://qed.econ.queensu.ca/pub/faculty/mackinnon/
The graphs are easy to interpret,
contain a considerable amount of information and even allow
direct comparison of statistics with different distributions under the
null.

Below is a simple program that demonstrates the use of the
main procedure, followed by the procedures themselves.  The code is
documented and hopefully user-friendly.  I'd welcome comments and
suggestions, or just information about similar procs that are already
out there.

Simon van Norden        Professeur invité
simon.van-norden@hec.ca or svn@alum.mit.edu or http://www.hec.ca/~p280
(514)340-6781 or fax (514)340-5632 or cell (514)993-0466
Service de l'enseignement de la finance, Ecole des H.E.C.
3000 Cote-Sainte-Catherine, Montreal QC, CANADA   H3T 2A7
***********************************************************/

new;
library pgraph;
graphset;

numtrial = 10000;
dta = rndn(numtrial,3);

dta[.,1] = dta[.,1].^2;
dta[.,2] = dta[.,2].^2 + dta[.,1];
dta[.,3] = dta[.,3].^2 + dta[.,2];

dof = {1 2 3};

_plegctl = 1;
fonts("simplex simgrma");

_plegstr =
"Ideal\000\202h\201[2](1)\000\202h\201[2](2)\000\202h\201[2](3)";
call pvalplot(dta,&pval_chi2,dof);

_plegstr =
"Ideal\000\202h\201[2](3)\000\202h\201[2](2)\000\202h\201[2](1)";
call pvalplot(dta[.,3 2 1],&pval_chi2,dof);

end;






/*        7 Feb 98 by Simon van Norden
        send comments & suggestions to svn@alum.mit.edu
PVALPLOT()

This program is designed for use with simulated data. Given a matrix of
data x, it generates
the data's theoretical p-values, G(x) (based on a user-provided
distribution.) It then constructs p-values based on the empirical CDF
for the data, F(x). Finally, it then generates two sets of graphs. The
first is a "P-Value Plot", which simply plots F(x) vs G(x).  If the
theoretical distribution is correct, we will see a 45-degree line. The
second is a "P-Value Discrepency Plot", which plots G(x)/F(x) vs G(x).
If the theoretical distribution is correct, we will see a horizontal
line at 1. Each type of graph is plotted for the full range of p-values
(0 to 1.0), and then for the range (0 to 0.1) most used for inference.

The reference for these graphs is Davidson and MacKinnon (1997).

The syntax of the program is
 pvalplot(x,&cdfproc,dof);
where
 x  is a matrix containing the statistics you want to examine
   (each row is one trial, each column is one statistic.)
 &cdfproc is a pointer to a procedure which calculates the *complement*
of the cdf
   (more on this below)
 dof  is a matrix containing the degrees of freedom parameters required
for
   &cdfproc

To use this procedure, you'll need a procedure to return the theoretical
CDF. This procedure will be called exactly once as
  dta = cdfproc(dta,dof);
Example of procs for the standard normal, the Chi-Squared and the F
distribution are included after the body of the main procedure. If you
need other distributions, it is up to you to provide a procedure.  Make
sure that
 a) it can evaluate all the columns and rows of dta simultaneously
 b) dof contains all the parameters which the cdfproc() might need


Before calling pvalplot, the user may want to use _plegctl and _plegstr
to set the legend for the graphs, or set _pframe = {0, 1}. The user may
control all aspects of the graphics aside from the following (which are
set by pvalplot()):
 ylabel
 xlabel
 window
 scale
 _ptitle
*/

proc pvalplot(dta,cdfproc,dof);

local cdfproc:proc, nobs, ncol, i, pval, xj;

@set up some constants@
nobs=rows(dta);
ncol=cols(dta);
pval=zeros(215,ncol);
xj = seqa(.001,.001,10)  | seqa(.015, .005,195) | seqa(.99, .001, 10);


@calculate p-values using theoretical distribution@
dta = cdfproc(dta,dof);

@calculate empirical CDF of p-values@
for i (1,ncol,1);
 pval[.,i] = counts(dta[.,i],xj);
endfor;
pval = cumsumc(pval)/nobs;

@Now do the graphics@

/* setup globals for both graphs */
xlabel("Theoretical p-value");

/* create 4 windows */
begwind;
window(2,2,0);

setwind(1);
_ptitle = "P-Value Plot";
scale(0|1, 0| 1);
ylabel("Empirical p-value");
xy(xj,xj~pval);

nextwind;
scale(0|0.1, 0| miss(1,1));
xy(xj[1:28,.],xj[1:28,.]~pval[1:28,.]);

_ptitle = "P Value Discrepancy Plot";
ylabel("Empirical / Theoretical p-value");

nextwind;
scale(0|1,0| miss(1,1));
xy(xj,ones(215,1)~(pval./xj));

nextwind;
scale(0|0.1,0| miss(1,1));
xy(xj[1:28,.],ones(28,1)~(pval[1:28,.]./xj[1:28,.]));

endwind;

retp(xj~pval);
endp;

proc pval_chi2(dta,dof);
 retp(cdfchic(dta,dof));
endp;

proc pval_normal(dta,dof);
 retp(cdfnc(dta));
endp;

proc pval_f(dta,dof);
 retp(cdffc(dta,dof[.,1],dof[.,2]));
endp;
