/* ==============================================================*/
/* The Ordinary Least Square (OLS) Module for GAUSS              */                      
/* ols.txt, Lastly Modified on March 25, 2001                    */
/* Programmed by Hun Myoung Park, kucc625@indiana.edu            */
/* Joint Ph.D. Student in Public Policy, Indiana University      */
/*                                                               */
/* Usage: {anova,beta,b_se,cov,corr,r,b_ci,cipi} = ols(x, y, t); */
/*        "#include ols.txt" should precede before being called  */    
/*                                                               */    
/* Input: x (n X (k-1) independent variable matrix)              */
/*        y (n X 1 dependent variable matrix)                    */
/*        t (t value for the degree of freedom, n-k, and alpha.  */
/*          if not specified, steps for ci and pi are ignored    */
/*                                                               */    
/* Output:anova (3 X 5 avova table matrix)                       */
/*               anova[2,3] -> "mse"; anova[2,4]->"see"          */
/*               anova[2,5] -> "R square";                       */
/*               anova[3,5] -> "Adj R square"                    */
/*        beta (k X 1 estimated coefficient, beta hat, matrix)   */
/*        b_se (standard error matrix for estimated coefficients */
/*        cov ( k X k covariance matrix)                         */
/*        corr (k X k correlation matrix)                        */
/*        r (n X 1 residual matrix)                              */
/*        b_ci (k X 5 confidence interval matrix                 */
/*                    for coefficients, betas)                   */
/*             b_ci[.,1]->t statistics; b_ci[.,2]-> p values     */ 
/*             b_ci[.,3]->lower ci; b_ci[.,4]->upper ci at alpha */ 
/*             b_ci[.,5]->difference between the lower and upper */ 
/*        cipi (n X 7 CI and predictioninverval matrix           */
/*                    for individual observations                */
/*             cipi[.,1]->y hat (estimated y values)             */ 
/*             cipi[.,2]->lower ci; cipi[.,3]->upper ci at alpha */ 
/*             cipi[.,4]->difference between the lower and upper */ 
/*             cipi[.,5]->lower pi; cipi[.,6]->upper pi at alpha */ 
/*             cipi[.,7]->difference between the lower and upper */ 
/* ==============================================================*/

PROC(8)=ols(x, y, t);

LOCAL n, k, h, df;
LOCAL vector, y_bar, v, beta, y_hat, r, anova;
LOCAL b_cov, b_se, b_t, b_corr, bc, h_diag, b_ci, n_cipi;

FORMAT/RD 15,5;

n=ROWS(x); h=1;        /* to calculte the number of observations */
vector=ONES(n, 1);                               /* n X 1 vector */
x=vector~x;                          /* merging (n X k) x matrix */
y_bar=meanc(y)*vector;                      /* n X 1 ybar matrix */

/*============= Calculating v, beta, y hat, residual =========== */
v=INV(x'*x);                          /* k X k covariance matrix */
beta=v*(x'*y);        /* Inv(X'X)*X'y to produce kX1 beta matrix */
y_hat=x*beta;            /* (n X k)*(k X 1) = n X 1 y hat matrix */
r=y-y_hat;                              /* n X 1 residual matrix */
k=ROWS(beta); df=n-k;        /* k=the total number of regressors */

/*================= Establishing ANOVA Table =================== */
anova=zeros(3,5);                    /* Establishing ANOVA Table */
anova[1,1]=k-h; anova[2,1]=df; anova[3,1]=n-h;            /* dfs */
anova[1,2]=(y_hat-y_bar)'*(y_hat-y_bar);                  /* ssr */
anova[1,3]=anova[1,2]/anova[1,1];              /* msr=ssr/(k-h); */
anova[2,2]=(y-y_hat)'*(y-y_hat);                          /* sse */
anova[2,3]= anova[2,2]/df;                     /* mse=sse/(n-k); */
anova[3,2]=(y-y_bar)'*(y-y_bar);                          /* sst */
anova[1,4]=anova[1,3]/anova[2,3];           /* f value= msr/mse; */
anova[1,5]=cdffc(anova[1,4], k-h, df);                /* p value */
anova[2,4]=SQRT(anova[2,3]);                  /* see = SQRT(mse) */
anova[2,5]=anova[1,2]/anova[3,2]; /* r square =ssr/sst=1-sse/sst */
anova[3,5]=1-((n-1)*(1-anova[2,5])*(1/df));  /*adjusted r square */

/* ======== Parameter Estimaters' Statistics ====================*/
b_cov=anova[2,3]*v;      /* s*s*INV(X'X) k X k covariance matrix */
b_se=sqrt(diag(v)*anova[2,3]);       /* standard error for betas */
b_t=beta./b_se;                        /* t value for H0: beta=0 */
bc=zeros(k, k);
b_corr=DIAGRV(bc, 1/b_se)*b_cov*DIAGRV(bc,1/b_se);/* correlation */   

/* ============ Confidence Interval Estimate ====================*/
b_ci=zeros(k, 5); n_cipi=zeros(n,7);

/* b_ci contains t, p, lower ci, upper ci, and ci diff. of beta  */
/* n_cipi contains y_hat, lower ci, upper ci, ci difference,     */ 
/*                 lower pi, upper pi, and pi difference         */

IF t > 0;                  /* Only when the t value is specified */

b_ci[.,1]=b_t;                                   /* t statistics */
b_ci[.,2]=2*cdftc(ABS(b_t), df);     /*1-PROBF(b_t#b_t, 1, df);  */
b_ci[.,3]=beta-t*b_se;                      /* lower ci of betas */
b_ci[.,4]=beta+t*b_se;                      /* upper ci of betas */
b_ci[.,5]=b_ci[.,4]-b_ci[.,3];         /* ci difference of betas */

/* 95% Confidence Intervals for the Mean of X */
h_diag=DIAG(x*v*x');   /* diagonal of "hat matrix"= X*INV(X'X)*X'*/
n_cipi[.,1]=y_hat;                               /* y hat matrix */
n_cipi[.,2]=y_hat-(t.*SQRT(h_diag*anova[2, 3]));      /*lower ci */
n_cipi[.,3]=y_hat+(t.*SQRT(h_diag*anova[2, 3]));      /*upper ci */
n_cipi[.,4]=2*(t.*SQRT(h_diag*anova[2, 3]));           /* 95% ci */

/* 95% Prediction intervals for each observation */
n_cipi[.,5]=y_hat-(t.*SQRT(h_diag*anova[2, 3]+anova[2, 3])); 
n_cipi[.,6]=y_hat+(t.*SQRT(h_diag*anova[2, 3]+anova[2, 3]));
n_cipi[.,7]=2*(t.*SQRT(h_diag*anova[2, 3]+anova[2, 3])); 

ENDIF;                                    /* End of IF statement */

/* ============ Returning the Vectors Calculated ================*/
RETP(anova, beta, b_se, b_cov, b_corr, r, b_ci, n_cipi);

ENDP;                                    /* The end of procedure */