* STABTEST.SRC.
*  Computes Hansen's stability test, from "Parameter Instability
*     in Linear Models", Journal of Policy Modeling, 14 (4), 1991, pp. 517-533.
*
* Procedure by Tom Maycock, Estima,
*   June 5, 1998
*
* Syntax:
*
*  @STABTEST([PRINT]/NOPRINT) depvar start end jointstat
*  # regressor list
*
* Computes individual and joint test statistics. Joint statistic
*  is stored in the "jointstat" parameter (if provided).
*
* If you want to be able to access the individual stats, just
*   declare the vector L as global, rather than local.
*
* The PRINT option controls whether or not the regression and
*  critical values table is printed.
*

procedure stabtest depvar start end jointstat
type series    depvar
type integer   start end
type real      *jointstat
option switch print 1

local series     e
local integer    m startl endl numobs entry var
local rect       xmat emat fmat smat jointv critval srow jointtemp
local vect[int]  varlist
local vect       s w l smat1
local real       jointsum

* Get regressor list:
enter(varying) varlist


dim critval(20,6)
compute critval = $
|| 0.748 , 0.593 , 0.470 , 0.398 , 0.353 , 0.240 | $
   1.07  , 0.898 , 0.749 , 0.670 , 0.610 , 0.469 | $
   1.35  , 1.16  , 1.01  , 0.913 , 0.846 , 0.678 | $
   1.60  , 1.39  , 1.24  , 1.14  , 1.07  , 0.883 | $
   1.88  , 1.63  , 1.47  , 1.36  , 1.28  , 1.08  | $
   2.12  , 1.89  , 1.68  , 1.58  , 1.49  , 1.28  | $
   2.35  , 2.10  , 1.90  , 1.78  , 1.69  , 1.46  | $
   2.59  , 2.33  , 2.11  , 1.99  , 1.89  , 1.66  | $
   2.82  , 2.55  , 2.32  , 2.19  , 2.10  , 1.83  | $
   3.05  , 2.76  , 2.54  , 2.40  , 2.29  , 2.01  | $
   3.27  , 2.99  , 2.75  , 2.60  , 2.49  , 2.22  | $
   3.51  , 3.18  , 2.96  , 2.81  , 2.69  , 2.41  | $
   3.69  , 3.39  , 3.15  , 3.00  , 2.89  , 2.59  | $
   3.90  , 3.60  , 3.34  , 3.19  , 3.08  , 2.77  | $
   4.07  , 3.81  , 3.54  , 3.38  , 3.26  , 2.95  | $
   4.30  , 4.01  , 3.75  , 3.58  , 3.46  , 3.14  | $
   4.51  , 4.21  , 3.95  , 3.77  , 3.64  , 3.32  | $
   4.73  , 4.40  , 4.14  , 3.96  , 3.83  , 3.50  | $
   4.92  , 4.60  , 4.33  , 4.16  , 4.03  , 3.69  | $
   5.13  , 4.79  , 4.52  , 4.36  , 4.22  , 3.86 ||

* Determine regression range:
inquire(reglist) startl>>start endl>>end
# varlist depvar

* Do regression:
linreg(print=print,def=stabeq) depvar startl endl e
# varlist
compute m = %nreg

* Note: To replicate Hansen's result, need e'e/T, rather than
*         e'e/(T-K), which is what RATS normally produces.
*       So, multiply %seesq by (T-K)/T to get Hansen results:
compute %seesq = %seesq*(%nobs-%nreg)/%nobs


* Store regressors into a rectangular array:
make xmat startl endl
# varlist

* Copy resids to an array:
make emat startl endl numobs
# e


* Compute V matrix
dim fmat(numobs,m+1)
do entry=1,numobs
 do i=1,m
  compute fmat(entry,i) = (xmat(entry,i)*emat(entry,1))
 end do i
 compute fmat(entry,m+1) = (emat(entry,1)**2 - %seesq)
end do entry
compute jointv = tr(fmat)*fmat

* Compute Joint stat
* First, build up S matrix for all entries (xe's in cols 1-m, e**2-s.e. in column m+1):
dim smat(entry,m+1)
do entry=1,numobs
 do var=1,m+1
   if entry==1
    compute smat(entry,var) = %if(var<=m, xmat(entry,var)*emat(entry,1), emat(entry,1)**2 - %seesq)
   else
    compute smat(entry,var) = %if(var<=m, smat(entry-1,var) + xmat(entry,var)*emat(entry,1), smat(entry-1,var) + emat(entry,1)**2 - %seesq)
 end do var
end do entry

* Now, sum S(t)*inv(V)*tr(S(t)) over all t:
dim srow(1,m+1)
compute jointsum = 0.0, jointtemp = 0.0
do entry=1,numobs
 ewise srow(i,j) = smat(entry,j)
 compute jointtemp = (srow*inv(jointv))*tr(srow)
 compute jointsum = jointsum + jointtemp(1,1)
end do

compute jointstat = jointsum/numobs

display 'Joint Statistic = '  ##.######## jointstat
display

* Compute Individual stability stats:
dim s(m+1)
compute s = %const(0.0)
dim w(m+1)
compute w = %const(0.0)
do entry=1,numobs
 do var=1,m+1
  compute s(var) = %if(var<=m, s(var) + xmat(entry,var)*emat(entry,1), s(var) + emat(entry,1)**2 - %seesq)
  compute w(var) = w(var) + s(var)**2
 end do
end do
dim l(m+1)
do var=1,m+1
 compute l(var) = (1.0/(numobs*jointv(var,var)))*w(var)
end

Display 'Individual Stability Test Statistics:'
compute eqnarray = %eqntable(stabeq)
do i=1,m
 if eqnarray(2,i)==0
   display @5 %label([series]eqnarray(1,i))  @30 ##.######## l(i)
 else
   display @5 %label([series]eqnarray(1,i)) @-1 '{' @-1 eqnarray(2,i) @-1 '}'  @30 ##.######## l(i)
end do
display @5 'Variance' @30 ##.######## l(m+1)

if PRINT
{
display
display 'Asymptotic Critical Values:'
display '---------------------------'
display '  Note: Use the values for one degree of freedom for individual statistics.'
display '        For example, reject stability hypothesis for an individual '
display '        component at 5% level if the statistic exceeds 0.470'
display
display @1 'Deg. of Freedom' @25 '1%' @35 '2.5%' @45 '5%' @55 '7.5%'  @65 '10%'  @75 '20%'
display @5 '1 '  @23 #.### critval(1,1) @34 #.###  critval(1,2) @43 #.### critval(1,3) @54 #.### critval(1,4) $
                     @63 #.### critval(1,5) @73 #.### critval(1,6)
* If degrees of freedom is > 1, print crit values for joint statistic:
if (m+1)>1
{
 display @5 (m+1) @23 #.### critval((m+1),1) @34 #.###  critval((m+1),2) @43 #.### critval((m+1),3) @54 #.### critval((m+1),4) $
                     @63 #.### critval((m+1),5) @73 #.### critval((m+1),6)
}
display
}

end procedure

