load auto1[134,8]=a:\auto1.asc;  /*  Loads an ascii matrix from the a: drive into gauss */
vnames=auto1[1,.];               /*  The first row has variable names in it */
                                 /*  Note, vnames contains words and is a string */ 
data=auto1[2:134,.];             /*  Rows 2 through 134 contain the numeric data */


timeline=seqa(1,1,133);   /*  A vector containing a sequence of observation numbers 1-133 */

price=data[.,6];      /*  Extracts the price series from the data matrix */
income=data[.,7];     /*  Extracts the income series from the data matrix */
gas=data[.,8];        /*  Extracts the gas series from the data matrix */

/* The quarterly dates and observation numbers are generated.
** This helps in matching dates with the correct observation
** number so that the correct subsamples can be generated.
*/  

d1=seqa(1959,1,34);
d1=d1.*.ones(4,1);
d2=seqa(1,1,4);
d2=ones(34,1).*.d2;
dates = 0 $+ ftocv(d1,2,0) $+ "." $+ ftocv(d2,1,0);

dates=dates[1:133,1];
xx=ftocv(seqa(1,1,133),1,0)~dates;
$xx;                 /* Prints observation numbers and dates to the screen */

/* Regression for sample 1959.1-1973.3 */

gas59_73=gas[1:59,1];
price59_73=price[1:59,1];

call ols(0,gas59_73,price59_73);

/* Regression for sample 1982.1 to 1992.1 */

gas82_92=gas[93:133,1];
price82_92=price[93:133,1];

call ols(0,gas82_92,price82_92);

/* Scatter plot corresponding to Figure 1.3.d in your book */
library pgraph;
_pstype=2;
xlabel("GAS");
ylabel("Price");
xy(gas82_92,price82_92);



/* Simple Least squares program */

y=gas82_92;                  /* Y is the new name for the dependent variable */
x1=price82_92;               /* X1 is the continuous independent variable    */
t=rows(y);                   /* Returns the number of observations           */
const=ones(t,1);             /* Generates a column of ones to estimate the intercept */

x=const~x1;                  /* Horizontally concatenates the ones and x1 to form X */
k=cols(x);                   /* Counts the number total number of independent variables */

b=invpd(x'*x)*x'*y;          /* Least Squares */
ehat=y-x*b;                  /* LS residuals  */
sse=ehat'*ehat;              /* Sum of squared errors */

sig2=sse/(t-k);              /*  Unbiased estimator of the variance of y */  

cov=sig2*invpd(x'*x);        /*  Estimates LS covariance matrix */
se=sqrt(diag(cov));          /*  Estimated standard errors */ 
tval=b./se;                  /*  T-ratios                  */ 

/*  Print out the estimates, std errors and t ratios */
xx=b~se~tval;
"      estimates         std errors        tvals ";
xx;