/*
** psn.sim
** (C) Copyright 1988-1993 by Aptech Systems, Inc.
** All Rights Reserved
**
** This command file simulates a data set with 3 independent
** variables and a dependent variable that is poisson distributed
** conditional on the x's.  The data set is saved as a GAUSS data set
** called psn.dat.
*/

    seed = 312453457;
    nobs = 100;
    b = { -1,
           1,
           1 };
    x = ones(nobs,1)~rndus(nobs,2,seed);
    yh = exp(x*b + rndns(nobs,1,seed));

/*
**  A random disturbance is added here even though that violates the
**  assumptions of the model that is estimated in max4.e and max5.e.
**  This is done here because in practice it is generally not possible to
**  to assume that all relevant exogenous variables are included in the
**  equation. Asymptotic theory shows that the estimates will be consistent,
**  and in addition that the correct standard errors may be computed from
**  the heteroskedastic-consistent covariance matrix of the parameters
**  (selected by setting __COVP = 3).
**
**  While optimizing the likelihood in max4.e or max5.e will
**  produce consistent estimates and standard errors when there is a
**  disturbance in the equation, more efficient estimates can be produced
**  by a correct specification of the log-likelihood.  This is done in
**  the COUNT program written by Gary King and available from APTECH SYSTEMS.
*/

    y = rndps(yh,seed);
    z = y~x;

    lbl = { Y, CNST, X1, X2 };

    create f0 = psn with ^lbl,0,8;
    if f0 == -1;                                                
        errorlog "Can't open output file";
        end;
    endif;
    if writer(f0,z) /= rows(z);
        errorlog "Disk Full";
        end;
    endif;
    f0 = close(f0);

/*
** This procedure will return a random poisson variable conditional
** on m.  The result will be conformable to m.  A missing value will
** be returned in any elements of m that are less than or equal to zero.
*/

proc rndps(m,seed);
    local s,n,r,c,i,j,w,k,z,f,m0,s0;

    r = rows(m);
    c = cols(m);
    s = rndus(r,c,seed);
    w = zeros(r,c);

    i = 0;
    do until i == r;
        i = i+1;
        j = 0;
        do until j == c;
            j = j+1;
            m0 = m[i,j];
            s0 = s[i,j];
            if m0 <= 0;
                w[i,j] = miss(0,0);
                continue;
            endif;
            k = 0;
            z = exp(-m0);
            f = z;
        A0:
            k = k+1;
            if z <= s0;
                f = f*(m0/k);
                z = z+f;
                goto A0;
            endif;
            w[i,j] = k;
        endo;
    endo;
    retp(w);
endp;