/* dec04.ex2                              */
 /*  analysis of California gas            */
 /*    consumption data                    */
 /*                                        */
  data a ;
    infile 'cagas.dat' ;
    input year qtr y ;
    q1 = (qtr = 1) ;                       * quarter dummies ;
    q2 = (qtr = 2) ;
    q3 = (qtr = 3) ;
    t+1 ;                                  * obs counter ;
	x0 = 1 ;
 run ;
  proc iml ;
    use a ;
    read all var{x0 t q1 q2 q3} into x ;
    read all var{y} into y ;
    print x y ;
    bign = nrow(y) ;                  * N obs ;
    p = ncol(x) ;                     * p cols ;
    bols = inv(x` * x) * (x` * y) ;   * usual bhat ;
    ehat = y - x * bols ;             * residuals ;
    sse = sum(ehat # ehat) ;
    elag = ehat[2:bign]// 0 ;
    rhohat = sum(elag # ehat) / sse ; * estimate lag correlation ;
    print rhohat ;
    v = do(1,bign,1) ;
    a = repeat(v,bign,1) ;
    b = abs( a - a`) ;                * (b)ij = abs(i-j) ;
    bigv = rhohat ## b ;
    xxi = inv(x` * x) ;
    covbols = (sse/bign)*xxi*(x`*bigv*x)*xxi ;
    print covbols ;
   quit ;
 run ;