Previous Page

|

Next Page

Example 8.3: Regression

This example shows a regression module that calculates statistics not calculated by the two previous examples. Here is the program.

  
       /*         Regression Routine               */ 
       /* Given X and Y, this fits Y = X B + E    */ 
       /* by least squares.                        */ 
  
    start reg; 
       n=nrow(x);                      /* number of observations */ 
       k=ncol(x);                         /* number of variables */ 
       xpx=x`*x;                               /* crossproducts */ 
       xpy=x`*y; 
       xpxi=inv(xpx);                   /* inverse crossproducts */ 
       b=xpxi*xpy;                        /* parameter estimates */ 
       yhat=x*b;                             /* predicted values */ 
       resid=y-yhat;                                /* residuals */ 
       sse=resid`*resid;                /* sum of squared errors */ 
       dfe=n-k;                      /* degrees of freedom error */ 
       mse=sse/dfe;                        /* mean squared error */ 
       rmse=sqrt(mse);                /* root mean squared error */ 
       covb=xpxi#mse;                 /* covariance of estimates */ 
       stdb=sqrt(vecdiag(covb));              /* standard errors */ 
       t=b/stdb;                        /* ttest for estimates=0 */ 
       probt=1-probf(t#t,1,dfe);     /* significance probability */ 
       print name b stdb t probt; 
       s=diag(1/stdb); 
       corrb=s*covb*s;               /* correlation of estimates */ 
       print ,"Covariance of Estimates", covb[r=name c=name] , 
              "Correlation of Estimates",corrb[r=name c=name] ; 
  
       if nrow(tval)=0 then return;   /* is a t value specified? */ 
       projx=x*xpxi*x`;                            /* hat matrix */ 
       vresid=(i(n)-projx)*mse;       /* covariance of residuals */ 
       vpred=projx#mse;       /* covariance of predicted values  */ 
       h=vecdiag(projx);                  /* hat leverage values */ 
       lowerm=yhat-tval#sqrt(h*mse); /* low. conf. lim. for mean */ 
       upperm=yhat+tval#sqrt(h*mse);     /* upper limit for mean */ 
       lower=yhat-tval#sqrt(h*mse+mse); /* lower limit for indiv */ 
       upper=yhat+tval#sqrt(h*mse+mse); /* upper limit for indiv */ 
       print ,,"Predicted Values, Residuals, and Limits" ,, 
       y yhat resid h lowerm upperm lower upper; 
    finish reg;
 

  
       /* Routine to test a linear combination of the estimates  */ 
       /* given L, this routine tests hypothesis that LB = 0.    */ 
  
    start test; 
       dfn=nrow(L); 
       Lb=L*b; 
       vLb=L*xpxi*L`; 
       q=Lb`*inv(vLb)*Lb /dfn; 
       f=q/mse; 
       prob=1-probf(f,dfn,dfe); 
       print ,f dfn dfe prob; 
    finish test; 
  
      /* Run it on population of U.S. for decades beginning 1790 */ 
  
    x= { 1 1 1, 
         1 2 4, 
         1 3 9, 
         1 4 16, 
         1 5 25, 
         1 6 36, 
         1 7 49, 
         1 8 64 }; 
  
    y= {3.929,5.308,7.239,9.638,12.866,17.069,23.191,31.443}; 
    name={"Intercept", "Decade", "Decade**2" }; 
    tval=2.57;  /* for 5 df at 0.025 level to get 95% conf. int. */ 
    reset fw=7; 
    run reg; 
    do; 
       print ,"TEST Coef for Linear"; 
       L={0 1 0 }; 
       run test; 
       print ,"TEST Coef for Linear,Quad"; 
       L={0 1 0,0 0 1}; 
       run test; 
       print ,"TEST Linear+Quad = 0"; 
       L={0 1 1 }; 
       run test; 
    end;
 

The results are shown in Output 8.3.1.