/******    Panel-Corrected Standard-Errors (Franzese,Mar96)   ****
"A GAUSS Procedure to Implement Beck-Katz PCSEs in Data Sets with
Non-Rectangular and/or Missing Data (with other Bells and Whistles)"

  It is a rare pooled data set whose data starts and ends in the same
time periods (t-p's) for each cross-section (c-s) and is devoid of missing
values.  Fortunately, this poses no intrinsic problem for the estimation of
the panel-corrected standard-errors (PCSEs) suggested by Beck and Katz
(APSR 1995).  Unfortunately, the RATS and GAUSS procedures written by
Beck and Katz to estimate PCSEs do assume rectangularity and the absence
of missing values.  In writing a GAUSS procedure to handle missing
values, I decided to automate the process as much as possible in the
hopes that this would save me some effort in the future (it already has)
and that it might be of some use to others working with this sort of data.
  The procedure takes as inputs vectors of variable and c-s names, the
dependent and independent variables, the number of c-s's, the first and
last t-p's, and scalar switches for optional (panel) WLS, and/or c-s and/or
t-p fixed effects.  Organize your data in the usual way (c-s 1, t-p 1 to
T; c-s 2, t-p 1 to T; etc).  Always include a constant as your first
independent variable and the dependent variable _name_ as your last
variable name.  (All this is noted in the procedure's comments.)  If the
data are non-rectangular, just put missing values in to rectangularize
it.  (Try to make your y-vector and x-matrix as compact as possible by
not including any t-p which no c-s uses; this will avoid procedure crashes.)
The output is simlar to the OLS.SRC (Aptech) output in GAUSS, except that
standard errors (and t-stats and p-levels) are from PCSEs and that the last
column (correlation of that X with Y) is replaced with the OLS standard errors.
  The only trick to dealing with missing data is to realize that the
cross-products of residuals, i.e., the sum over t of e(i,t)*e(j,t) does
not contain the same number of valid observations for each i and j.  Thus
the denominator of each element of the variance-covariance matrices are
(potentially) different.  The procedure offered here handles this by
creating a vector of ones for valid and zeros for missing observations.
This vector can be reshaped and/or manipulated (summed, multiplied, ...) as
necessary to divide each element of the variance-covariance matrices by
the correct number of valid observations.  (The code is heavily remarked
and should be understandable if you know GAUSS and PCSEs.)
  There are three additional bivariate options I have included in the
procedure.  First, you can choose WLS or OLS.  If WLS is requested, OLS
is run first (this first-stage regression is currently set not to be
output to the screen; you may change it by hardwiring the appropriate
_output=0 line to _output=1).  The residuals from that are squared and
regressed on c-s dummies (in such a way that the F-stat for the
regression is a test of the panel-weight model against homoskedasticity).
(This regression is currently set to be output to the screen; you may
change that by hardwiring the appropriate line.)  The inverse of the
square roots of the fitted values from that regression are your usual
panel weights.  (An alternative model of residual variance can be
substituted by hardwiring changes to the lines noted in the procedure.)
WLS is then performed by running OLS on the transformed data, and PCSEs
are calculated from that weighted regression.  Weighted and unweighted
statistics are reported.
  The procedure will also automatically create a set of c-s and/or t-p
dummies for fixed effect models if requested to do so.  For a recent
discussion of when this is appropriate and of interpretation see Smith
(1995).  (Keep the constant in your independent variables; it'll be
deleted for you.  Make sure you do not otherwise create perfect
colinearity or the procedure with end with an error message from OLS.)
  Finally, in writing this procedure, I noticed that OLS.SRC ((c)
APTECH), generates incorrect DW statistics in the presence of missing
values (it "packs" vectors with missing values, making some of the
"adjacent" observations not truly temporally adjacent; APTECH has
been notified of the error).  This has been corrected in my procedure,
but be warned that whenever you are using the DW statistic produced by
GAUSS OLS.SRC outside of PCSE.G you are not getting the correct DW.
(Incidentally, the DW is invalid when there are lagged dependent
variable(s) in the equation.  An alternative test for
first-order serial correlation with a lag is:

h = (approx) (1-DW/2)*[n/(1-n*V(b1))] which is distributed
assymptotically standard normal (where V(b1) is the estimated variance of
the coefficient on the lagged dependent variable).

Q and LM tests are other valid alternatives, see Johnston 1984 sections
8-5 through 10-3.)

  Any suggestions on improving this procedure would be gratefully
received.  The procedure is public domain and you may use it as you wish
(though a citation of the source would be appropriate/appreciated).  I am
confident the procedure correctly calculates PCSEs; however, I cannot
accept responsibility for any difficiulties, financial or otherwise,
arising from its use.  (My lawyer/sister insisted I officially note that.)

*****************************************************************

A debt of gratitude is owed Curtis Signorino (Harvard University) for his
assistance in optimizing the procedure for memory and time use.

References:

Beck, N. and J. Katz, "What to do (and not to do) with
Time-Series-Cross-Section Data in Comparative Politics," APSR v. 89 n. 3,
September 1995, pp. 634-47.

Johnston, J., _Econometric Methods_, 3rd ed., New York: McGraw Hill, 1984.

Smith, M., "Time-Serial Designs and Cross-Sectional Designs: Uncovering
the Structural Logic of Pooled Analyses," Paper presented at the 1995 APSA.
***********************************************************************/
/***************  PROGRAM CODE FOLLOWS **************************/
@ Gauss Procedure to Estimate Linear Regression, Implementing
* Panel-Corrected Standard-Errors (PCSE) (Beck and Katz (APSR 1995))
* in the Presence of Non-rectangular and/or Missing Data
*
* Please cite "PCSE.G: A Gauss Procedure to Implement Panel-Corrected
* Standard-Errors in Non-Rectangular Data Sets."
*
* by Robert J. Franzese, Jr. (Harvard University) October 1995
* Revised 3/19/96
*
* Code is written and submitted for public, non-commercial use.
* The authors assumes no responsibility whatsoever for the
* performance of the procedure.
*
* Procedure handles non-rectangular data and/or missing observations,
* and (optional) creates and uses cross-section and/or time-period
* fixed effects, and (optional) estimates and uses panel weights (other
* linear variance models can be estimated if desired by editing (hardwiring)
* the line `vx = ...' (remarked in the procedure))
*
* AUTHOR'S NOTE:
* The assistance of Curt Signorino (Harvard University) in optimizing
* the procedure to reduce computer time and memory use is gratefully
* acknowledged.  The author was funded by a Mellon Foundation Dissertation
* Completion Grant during the writing of the procedure.
*
* FORMAT:
* {m,b,stb,zz,stderr,se,sigma,cx,rsq,rbsq,resid,dwstat,nobs,degfree,
*  tstats,poft,ursq,urbsq,usigma,udw,ures,corre}
* =pcse(vnames,cnames,y,x,nocs,tstart,tend,weight,cdum,tdum);
*
* Procedure assumes data organized {[c1,t1],[c1,t2],...,c[2,t1],...,etc.}',
* that each cross-section has the same number of rows of the data matrix
* allotted to it, though some of those rows may be missing data, and that
* the first column of "indvars" is the constant.  (Thus, if your data is non-
* rectangular, you must rectangularize it by putting missing values in. Any
* observations you wish to drop must be fed as missing in the dep var or at
* least one indep var.)
*
* Procedure returns: 22 total, in the following order:
    * m = (matrix) moment matrix (OLS/WLS)
    * b = (vector) coefficients (OLS/WLS)
    * stb = (vector) standardized coefficients (OLS/WLS)
    * zz = (matrix) variance-covariance of coefficients (OLS/WLS + PCSE)
    * stderr = (vector) coefficient standard-errors (OLS/WLS)
    * se = (vector) coefficient standard-errors (OLS/WLS + PCSE)
    * sigma = (scalar) std err of the estimate (WLS, same as usigma if OLS)
    * cx = (matrix) variance-covariance of X and Y (OLS/WLS)
    * rsq = (scalar) R-squared (WLS, same as ursq if OLS)
    * rbsq = (scalar) R-bar-squared (WLS, same as urbsq if OLS)
    * resid = (vector) residuals (WLS, same as ures if OLS)
    * dwstat = (scalar) Durbin-Watson statistic (WLS, same as udw if OLS)
    * nobs = (scalar) number of valid observations
    * degfree = (scalar) degrees of freedom of coefficient estimates
    * tstats = (vector) t-stats of test that coeffs=0 (OLS/WLS + PCSE)
    * poft = (vector) 2-sided prob of those t-tests (OLS/WLS + PCSE)
    * ursq = (scalar) R-squared (unweighted residuals used)
    * urbsq = (scalar) R-bar-squared (unweighted residuals used)
    * usigma = (scalar) std error of estimate (unweighted residuals used)
    * udw = (scalar) Durbin-Watson statistic (unweighted residuals used)
    * ures = (vector) unweighted residuals
    * corre = (matrix) Cross-section Correlation of Residuals

* Procedure Inputs: 10 in the following order:
    * vnames = character vector: variable names, last being dep var
    * cnames = character vector: cross-section unit names
    * y = vector: observations on dependent variable
    * x = matrix: obs. on independent variables including constant (first)
        * time and or c-s dummies can be added for you, omit them here
    * nocs = scalar: number of cross-sectional units
    * tstart = scalar: first time period (earliest from all c-s units)
    * tend = scalar: last time period (latest from all c-s units)
    * weight = scalar: 1 for WLS, 0 for OLS
        * default WLS model is panel weights, can change that directly by
        * hard-wiring the line `vx = ...' in this procedure
    * cdum = scalar: 1 = add c-s dummies, 0 = no c-s dummies
    * tdum = scalar: 1 = add time dummies, 0 = no yime dummies
        * if tdum=1 and cdum=1, procedure drops 1st time dummy @

proc (22) =
    pcse(vnames,cnames,y,x,nocs,tstart,tend,weight,cdum,tdum);

local   i,j,t,cdums,tdums,denom,vcve,corre,zz,cm,valid,nindvar,nobs,
        wx,wy,vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat,se,degfree,
        tstats,poft,results,mask,fmt,ures,ursq,urbsq,usigma,udw,tv,rbsq,
        sse,usse,nots,res2,vx,vare,minvar,wvare,wtv,res,resm,uresm,xm,ym,
        wxm,wym,tseq,rresm,rvalid,diagvcve,uresm1l,validm1l,resm1l;

/*declare vars local to procedure*/

format 12,6; /* resets default format for screen output */

nots = tend - tstart + 1; /* computes number of time periods */

if rows(y) /= rows(x); /* this set of ifs checks if input valid */
    "ERROR: Different number of rows in X and Y";
    end;
elseif rows(y) /= nocs*nots;
    "ERROR: Rows in Y and X do not equal # of cs times # of t-p";
    end;
elseif rows(cnames) /= nocs;
    "ERROR: Number of c-s names does not equal number of c-s";
    end;
elseif rows(vnames) /= cols(x)+1;
    "ERROR: Number of variable names does not equal columns X + columns Y";
    end;
else; /* no errors in input detected, create an indicator for valid obs */
    valid = y~x; /*make valid the data matrix */
    valid = (valid./=miss(1,1)); /*set each element to 1 if not miss 0 else*/
    valid = prodc(valid'); /*creates column vector with 1 if valid obs */
    validm1l = missrv(lag1(valid),0); /* lag of valid */
    y = valid.*y+(1-valid).*"miss"; /*put "miss" in y for all missing*/
    y = miss(y,"miss") ; /* convert "miss" to missing value code */
    x = valid.*x+(1-valid).*("miss".*ones(1,cols(x))); /* ditto for x */
    x = miss(x,"miss") ;
endif;

/** Describing of Data Set Shape **/
nobs=sumc(valid); /* compute number of valid observations */
print ftos(nots,"Maximum Number of Time-Periods   = %*.*lf",0,0);
print ftos(nocs,"Maximum Number of Cross-Sections = %*.*lf",0,0);
print ftos(nots*nocs,"NoT x NoC (possible observations) = %*.*lf",0,0);
print ftos(nobs,"Actual Number of Observations Valid  = %*.*lf",0,0);
print ftos(nots*nocs-nobs,"Number of Invalid Observations   = %*.*lf",0,0);

if cdum==1; /* create a matrix of c-s dummies if requested */
            /* delete constant and adjust variable names to match */
    vnames = vnames[2:rows(vnames),.]; /* delete var name for constant*/
    vnames = cnames|vnames; /* append c-s names to remaining var names */
    cdums = (eye(nocs).*.ones(nots,1)); /* create matrix of c-s dummies */
    cdums = valid.*cdums+(1-valid).*("miss".*ones(1,nocs)); /*put miss in*/
    cdums = miss(cdums,"miss");
    x=x[.,2:cols(x)]; /* remove constant from x matrix */
    x=cdums~x; /* c-s dummies are put at front of x-matrix */
elseif cdum /= 0;
    "ERROR: `cdum' argument must be a scalar, 0 (off) or 1 (on)";
    end;
endif;

/** generate time-period dummies **/

if tdum==1; /* create a matrix of time dummies if requested */
            /* delete constant if cdum hasn't already done so */
            /* append time dummies matrix to front of x and */
            /* adjust vnames to match the resulting x matrix */
    tdums = (ones(nocs,1).*.eye(nots)); /* create matrix of t-p dummies */
    tdums = valid.*tdums+(1-valid).*("miss".*ones(1,nots)); /* put miss
in */
    tdums = miss(tdums,"miss");
    tseq = seqa(tstart,1,(tend-tstart+1));
    if cdum == 0 ;
        vnames = vnames[2:rows(vnames),.]; /* remove constant name */
        x = x[.,2:cols(x)]; /* remove constant */
        vnames = (0$+"T"$+ftocv(tseq,0,0))|vnames; /*append tnames to
v-names*/
        x = tdums~x; /* append tdums to x-matrix */
    else ; /* if c-s dums, same, except trim first time-period dummy */
        vnames = trimr((0$+"T"$+ftocv(tseq,0,0)),1,0)|vnames;
        x = (trimr(tdums',1,0))'~x;
    endif;
elseif tdum /=0;
    "ERROR: `tdum' argument must be a scalar, 0 (off) or 1 (on)";
    end;
endif;

nindvar=cols(x); /* now compute total number of independent variables */
print ftos(nindvar,"Number of Independent Variables = %*.*lf",0,0);
"";
degfree = nobs-nindvar; /* compute degrees of freedom */

if weight==0; /* if weighting not requested, do OLS */
    "Using OLS in First and Second Stage...";
    "";
    __altnam = vnames; /* feeds vnames to OLS */
    __miss = 1; /* instructs OLS to omit missing data */
    __con = 0; /* instructs OLS not to add another constant */
    _olsres = 1; /* instructs OLS to compute residuals and DW-stat */
    __output = 0; /* set this to 1 if you want OLS results printed out */

{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat }=OLS(0,y,x); /* OLS */

    rbsq = 1-(1-rsq)*((nobs-1)/degfree); /* compute r-bar-squared */
    res = y-x*b; /* compute nocs*nots by 1 vector of residuals */
    resm = missrv(res,0); /* an alternate resid vector with 0 for missing*/
    resm1l = missrv(lag1(resm),0); /* lag of resm */
    ures=res; /* unweighted residuals, same as OLS here */
    uresm = missrv(ures,0); /* alternate with zeros for missings */
    sse=resm'resm; /* compute sum of squared errors */
    usse=sse; /* unweighted sum of squared errors, same as OLS here */
    ursq=rsq; /* unweighted r-squared, same as OLS here */
    urbsq=rbsq; /* unweighted r-bar-squared, same as OLS here */
    usigma=sigma; /* unweighted std err of est, same as OLS here */

/** nb: GAUSS OLS computes incorrect DW when missing values present **/

    denom = sumc(valid.*validm1l.*resm^2); /* denom for DW is sum sq res */
    dwstat = sumc(valid.*validm1l.*(resm-resm1l)^2); /* numerator is... */
    dwstat = dwstat/denom; /*...sum of squared differences */
    udw=dwstat; /* unweighted DW-stat, same as OLS here */

    "First-stage done, now computing Panel V-Cov of Residuals...";

         @ covariance of resid[i,t] with resid[j,t] is computed
         * taking care to divide sum of cross-products only by the
         * number of valid observations @

    rvalid = reshape(valid,nocs,nots)'; /* reshape valid to nots x nocs */
    rresm = reshape(resm,nocs,nots)'; /* ditto for residual vector */
    vcve = rresm'*rresm; /* generate nocs by nocs matrix of crossprods */
    denom = rvalid'*rvalid;
    denom = miss(denom,0);
    if ismiss(denom);
         "ERROR: A c-s unit exists w/o any observations or w/o "\
         "any observations in common with another c-s unit.";
         "You must remove that c-s unit from the data passed to PCSE.";
         end;
    endif;
    vcve = vcve./denom;

    "Done, now computing cross-section residual-correlation matrix...";

    diagvcve = sqrt(diag(vcve));
    corre = vcve./(diagvcve*diagvcve');

    "Done, now computing PCSE V-COV of b...";
    xm = missrv(x,0); /* GAUSS doesn't like to do m-algebra with missing */
    ym = missrv(y,0); /* obs, oblige it by replacing with 0.  This has */
                      /* no effect on x'x, just watch the denominators */
    tv = ym'ym-nobs*((sumc(ym)/nobs)^2); /* sum of squared variance of y */
    wtv = tv; /* sum of squared variance of weighted y, same for OLS */

    zz = (inv(xm'xm))*xm'(vcve.*.(eye(nots)))*xm*(inv(xm'xm));
    /*PCV-COV(b)*/

    se = sqrt(diag(zz)); /* PCSE's */
    tstats = b./se; /* PC t-stats */
    poft = 2*cdftc(abs(tstats),degfree); /* 2-sided p of |PCSE t's| */

elseif weight==1; /* if weighted least squares requested */
    "Using OLS in First Stage...";
    "";
    __altnam = vnames; /* feed OLS the relevant variable names */
    __miss = 1; /* instruct OLS to omit missings */
    __con = 0; /* instruct OLS not to add another constant */
    _olsres = 1; /* instruct OLS to compute residuals and DW-stat */
    __output = 0; /* change this to 1 if you want the OLS output */

{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat }=OLS(0,y,x); /* OLS */

"Done, now calculating weighting series using OLS on OLS residuals squared";

    res2=(y-x*b).*(y-x*b); /* compute nocs*nots by 1 vector of residuals */
    if cdum==0; /* check to see if matrix of c-s dummies exists, if not...*/
        cdums = (eye(nocs).*.ones(nots,1)); /* create matrix of c-s dums */
        cdums = valid.*cdums+(1-valid).*("miss".*ones(1,nocs));
        cdums = miss(cdums,"miss");
    endif;

@ Residual Variance Model can be changed by editing the line `vx = ...' @

    vx = ones(nocs*nots,1)~cdums[.,2:nocs] ; /* ind vars in res var model */

@ The following process creates a weighting series by regressing the squares
* of the OLS residuals on the cross-section dummies (and whatever other
* explanitory variables you might have included in the preceding line.)
* The inverse of the square roots of the fitted values are the weights.
* Without other explanitors than the c-s dummies, this is exactly equivalent
* to panel-weighted least-squares @

    "----> n.b., F-stat is test of variance model against homoskedasticity";
    "";

@ if `vx= ...' is changed, the following line must be altered to match @

    __altnam = "one"|cnames[2:nocs,1]|"SqOLSres";
    __output = 1; /* 0 if you don't want the res var model output */

    { vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat }=OLS(0,res2,vx);

    vare = vx*b; /* make nocs*nots by 1 vector of fitted res variance */

@ The following if...endif checks for negative fitted values
* in the variance model, replacing them with the minimum positive value.
* (author's note: substitution of a more appropriate model than a linear
* model of variance may be particularly appropriate if negative fitted
* values are common) They are not possible in the panel-weight set-up. @

"";
"Calculating Weighting Series..."; /* Set negative wts to min positive */
    if minc(vare) <= 0;
        minvar = minc((vare.>0).*vare+(vare.<=0).*maxc(vare));
        vare = (vare.>0).*vare+(vare.<=0).*minvar;
    endif;

    wvare = ones(rows(vare),1)./sqrt(vare); /* weighting series */
    wx = wvare.*x; /* transform x to weighted x */
    wy = wvare.*y; /* transform y to weighted y */
    wxm = missrv(wx,0); /* create alternates wx and wy where miss are 0 */
    wym = missrv(wy,0);

"Done, now doing WLS using OLS on transformed variables...";

    __altnam = vnames;
    __output = 0; /* change this to 1 if you want WLS output */

{ vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat }=OLS(0,wy,wx); /*WLS*/

    res = wy-wx*b; /* nocs*nots vector of weighted residuals */
    resm = missrv(res,0); /* alternate with 0 for missing */
    resm1l = missrv(lag1(resm),0); /* lag of resm */
    rbsq = 1-(1-rsq)*((nobs-1)/degfree); /* weighted r-bar-squared */
    ures = y-x*b; /* calculate unweighted residuals */
    uresm = missrv(ures,0); /* alternate, 0 for missings */
    ym = missrv(y,0); /* alternate y with zeros for missings */
    tv = ym'ym-nobs*((sumc(ym)/nobs)^2); /* unweighted sum of var(y) */
    wtv = wym'wym-nobs*((sumc(wym)/nobs)^2); /* weighted sum of var(y) */
    ursq = (tv - uresm'uresm)/tv; /* unweighted r-squared */
    urbsq = 1-(1-ursq)*((nobs-1)/degfree); /* unweighted r-bar-squared */
    sse=resm'resm; /* sum of wwighted residuals variance */
    usse = uresm'uresm; /* unweighted sum of squared residuals */
                 /* n.b.,in WLS, unweighted resids may not have mean 0 */
    usigma = sqrt((uresm'uresm)/degfree); /* unwtd std err of estimate */
            /* computing unweighted DW-stat */
    uresm1l = lag1(uresm); /* create lagged residual (alternate) */
    uresm1l = missrv(uresm1l,0); /* replace missings with zeros */
    validm1l = missrv(lag1(valid),0); /* ditto for valid */
    denom = sumc(valid.*validm1l.*uresm^2); /* denom for DW is sum sq res */
    udw = sumc(valid.*validm1l.*(uresm-uresm1l)^2); /* numerator is... */
    udw = udw/denom; /*...sum of squared differences */

/** nb: GAUSS OLS computes incorrect DW when missing values present **/

    denom = sumc(valid.*validm1l.*resm^2); /* denom for DW is sum sq res */
    dwstat = sumc(valid.*validm1l.*(resm-resm1l)^2); /* numerator is... */
    dwstat = dwstat/denom; /*...sum of squared differences */

"Done, computing PCSE V-COV of WLS residuals...";
    rvalid = reshape(valid,nocs,nots)'; /* reshape valid to nots by nocs */
    rresm = reshape(resm,nocs,nots)'; /* ditto for resids */
    vcve = rresm'*rresm; /* calculate nocs by nocs cross-prod of resids */
    denom = rvalid'*rvalid; /* calculate nobs for each element of vcve */
    denom = miss(denom,0);
    if ismiss(denom);
         "ERROR: A c-s unit exists w/o any observations or w/o "\
         "any observations in common with another c-s unit.";
         "You must remove that c-s unit from the data passed to PCSE.";
         end;
    endif;
    vcve = vcve./denom; /* cross-prods over nobs for each is covar matrix */

"Done, computing cross-section residual-correlation matrix...";

    diagvcve = sqrt(diag(vcve)); /* sq root of diag elements = std dev */
    corre = vcve./(diagvcve*diagvcve'); /* vcv over sd()sd() = corr */

"Done, computing PCSE V-COV of b..."; /* PCSE V-COV(b) with wtd resids */

    zz = (inv(wxm'*wxm))*wxm'*(vcve.*.(eye(nots)))*wxm*(inv(wxm'*wxm));
            /* PCV-COV */

    se = sqrt(diag(zz)); /* PCSE's from WLS */
    tstats=b./se; /* PC t-stats from WLS */
    poft = 2*cdftc(abs(tstats),degfree); /* 2-sided p of |PC t| */

else;

    "Weight parameter must be set as 0 (=off) or 1 (=on)";
    end;

endif;
"** PCSE Results **"; /* printing out the results pretty-like */
"";
print ftos(nobs,"Valid Cases:  %*.*lf",20,0);;
print ftos(vnames[rows(vnames)],"      Dependent Variable:%*.*s",20,0);
print ftos(nocs*nots-nobs,"Missing Cases:%*.*lf",20,0);;
print "      Deletion Method:               ";;
print "Listwise"; /** always uses listwise **/
if weight == 1;
    print ftos(wtv,"Wtd Total SS: %*.*lf",20,3);;
    print ftos(degfree,"      Degrees of Freedom:%*.*lf",20,0);
endif;
print ftos(tv,"Total SS:     %*.*lf",20,3);;
print ftos(degfree,"      Degrees of Freedom:%*.*lf",20,0);
if weight == 1;
    print ftos(rsq,"Wtd R-Squared:%*.*lf",20,3);;
    print ftos(rbsq,"      Wtd Rbar-Squared:  %*.*lf",20,3);
endif;
print ftos(ursq,"R-Squared:    %*.*lf",20,3);;
print ftos(urbsq,"      Rbar-Squared:      %*.*lf",20,3);
if weight == 1;
    print ftos(sse,"W Residual SS:%*.*lf",20,3);;
    print ftos(sigma,"      Wtd Std Error Est: %*.*lf",20,3);
endif;
print ftos(usse,"Residual SS:  %*.*lf",20,3);;
print ftos(usigma,"      Std Error of Est:  %*.*lf",20,3);
print ftos(udw,"DW-Stat:      %*.*lf",20,3);;
print ftos(dwstat,"      Wtd DW-Stat:       %*.*lf",20,3);
print;
print "                         PCSE Std                 Prob    Sta"\
"ndardized  Reg Std";
print "Variable     Estimate     Error      t-value      >|t|      E"\
"stimate      Error";
print "-------------------------------------------------------------"\
"------------------";
results = vnames[1:nindvar,1]~b~se~tstats~poft~stb~stderr;
mask = 0~1~1~1~1~1~1;
let fmt[7,3] = "-*.*s" 9 8 "*.*lf" 12 6 "*.*lf" 12 6 "*.*lf" 12 6 ""\
        "*.*lf" 10 4 "*.*lf" 12 6 "*.*lf" 12 6;
call printfm(results,mask,fmt);
"";

retp(m,b,stb,zz,stderr,se,sigma,cx,rsq,rbsq,res,dwstat,
        nobs,degfree,tstats,poft,ursq,urbsq,usigma,udw,ures,corre);

endp;

