/* ==============================================================*/
/* The White Heteroscedasticity Consistent Estimator Module      */
/* wgls.txt (GAUSS), Lastly Modified on March 25, 2001           */
/* Programmed by Hun Myoung Park, kucc625@indiana.edu            */
/*                                                               */
/* Usage : {w_cov, w_se} = wgls(x, beta, r);                     */
/*        "#include wgls.txt" should precede before being called */    
/*                                                               */
/* Input : x ( n X (k-1) independent variable matrix)            */
/*         beta (k X 1 estimated coefficient matrix)             */
/*         r (n X 1 residual matrix)                             */
/*                                                               */
/* Output: w_cov is the k X k covariance matrix                  */
/*         w_se is the matrix for se, t, p values                */
/*         w_se[.,1] Standard errors of the coefficients         */
/*         w_se[.,2] T statistics of the coefficients (for H0)   */
/*         w_se[.,3] P values                                    */
/* ==============================================================*/

PROC(2) = wgls(x, beta, r);

LOCAL n, k, i;                                        /* scalars */
LOCAL vector, w_cov, w_se;                   /* matrix variables */
FORMAT/RD 15,5;                             /* to specify format */

n=ROWs(x); vector=ONES(n, 1);
x=vector~x; k=COLS(x);

             /* calculate estimated asymptotic covariance matrix */
i=1; w_cov=ZEROS(k, k);
DO WHILE i<=n;
    w_cov= w_cov + (r[i]'*r[i])*(x[i,.]'*x[i,.]);
    i=i+1;
ENDO; 
                              /* estimated asymptotic covariance */
w_cov=n*INV((x'*x))*(w_cov/n)*INV((x'*x));  

w_se=ZEROS(k,3);           /* to store standard errors, t, and p */
w_se[.,1]=SQRT(DIAG(w_cov));     /* to calculate standard errors */
w_se[.,2]=beta./w_se[.,1];          /* to calculate t statistics */
w_se[.,3]=2*CDFTC(ABS(w_se[.,2]), n-k);              /* p values */

RETP(w_cov, w_se);

ENDP;