/* ==============================================================*/ /* 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 */