/************************************************
********	240C - GAUSS - Example 1	*********
*************************************************
This program generates exponential random
variables and checks that the LLN and the CLT
apply for increasingly larger samples
************************************************/
new;
ssize   = 2|10|50|100|500;      @ ssize is a vector of sample sizes@
nss     = rows(ssize);          @ number of smaple sizes @
mcreps  = 5000;                 @ Monte Carlo replications @

lambda  = 5;                    @ Mean of the exponential r.v. @

/************************************************************************
Generate exponential random variables
*************************************************************************
General Comments:

1.- GAUSS does not have an exponential r.v. generator. Hence, one has to
generate uniform r.v.s and then do an appropriate transformation

2.- Notice that we will have to do a loop to generate samples of size
equal to each of the elements of the vector ssize and repeated mcpreps.

3.- Any matrices on which you want to store information sequentially
element by element, require that you initialize them so that GAUSS knows
where to input your entries.

4.- The following loop will calculate the sample mean per MC replication,
and then the SD, skewness, kurtosis, and Jarque-Bera statistic of the 
sample mean from the Monte Carlo replications. Note that in practice we
only observe one sample. I do this to get at the distribution artificially
for the example
***************************************************************************
Initialize the matrices that store the desired information
**************************************************************************/

mx  =   zeros(nss,1);     @ stores the sample mean for each sample size in ssize @
sx  =   zeros(nss,1);     @ stores the standard deviation of the sample mean across Monte Carlo replications @
hmx =   zeros(mcreps,nss);@ stores sample means to do graphs @

/*************************************************************************
Set up the loop
*************************************************************************/
i = 1;                  @ iteration counter @

do while i le nss;

    u = rndu(ssize[i],mcreps);      @ generate the uniform random variables @
    x = -lambda*ln(u);              @ transform from uniform into exponential r.v. @
    
    hmx[.,i]=   meanc(x);           @ sample mean for each of the mcreps @
    mx[i]   =   meanc(meanc(x));    @ take the mean of the mcreps sample means @
    sx[i]   =   stdc(meanc(x));     @ take the SD of the mcreps sample means @

    i = i+1;
endo;

sers = "Average Sample Mean "$~"Standard Deviation ";

/***************************************************************************
Output to screen
***************************************************************************/
sers;
mx~sx;

gtitles = "RAW DATA"$~"Sample Size = 2"$~"Sample Size = 10"$~"Sample Size = 100"$~"Sample Size = 1000"$~"Sample Size = 5000";
begwind;
_ptitlht = 0.25;
window(2,3,0);

setwind(1);
title(gtitles[1]);
{b,m,f} = hist(x,50);

i = 1;
do while i le nss;

	nextwind;
	title(gtitles[i+1]);
	{b,m,f} = hist(hmx[.,i],25);
	i = i+1;
endo;
endwind;

// The following command saves your graph to a postcript file
ret=tkf2ps("graphic.tkf","graphic.ps");


/*******************************************************************
load the libraries
*******************************************************************/
library pgraph;
graphset;