*
* STABTEST.SRC.
*  Computes Hansen's stability test, from "Parameter Instability
*     in Linear Models", Journal of Policy Modeling, 14 (4), 1991, pp. 517-533.
*
* Syntax:
*
*  @STABTEST( options ) depvar start end
*  # regressor list
*
* Computes individual and joint test statistics.
*
* Options:
*   [PRINT]/NOPRINT
*     Controls whether or not the regression and stability test results
*     are printed
*   SMPL=SMPL Series
*
* Revision Schedule:
*   06/1998 Written by Tom Maycock, Estima
*   08/2003 Revised by Tom Doan, Estima to handle missing values,
*      interpolate p-values.
*
function HTestCrit test m
type real HTestCrit
type real test
type integer m
*
* HTestCrit returns an approximate p-value for Hansen's
* stability statistic. For numbers between the tabled values
* it does a linear interpolation. It returns 0 for those below
* the bottom significance level (.01) and 1 for those above
* the top significance level (.20).
*
local rect critval
local vect pvalues
local integer i
local real p
local string htemp
*
dim critval(20,6) pvalues(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 ||
compute pvalues=||.010,.025,.050,.075,.100,.200||
*
compute p=1.0
do i=1,6
   if test>critval(m,i) {
      if i==1 {
         compute p=0.0
         break
      }
      else {
         compute p=pvalues(i-1)+(pvalues(i)-pvalues(i-1))*(test-critval(m,i-1))/(critval(m,i)-critval(m,i-1))
         break
      }
   }
end do i
compute HTestCrit=p
end
*
procedure stabtest depvar start end
type series    depvar
type integer   start end
*
option switch print 1
option series smpl

local series     e smpll
local integer    i startl endl time
local vector     f s x
local symm       v vinv
local real       sigsq
*
declare vector %HSTATS
declare real   %HJOINT
*
* Run the regression. In the process, this gets the supplementary
* card for the regressors.
*
linreg(print=print,smpl=smpl,define=stabeq) depvar start end e
compute startl=%regstart(),endl=%regend()
set smpll startl endl = %valid(e)
*
* Hansen uses e'e/T for the estimate of the variance
*
compute sigsq = %rss/%nobs
*
dim f(%nreg+1) s(%nreg+1) v(%nreg+1,%nreg+1)
compute v=%const(0.0)
do time=startl,endl
   if .not.smpll(time)
      next
   compute x=%eqnxvector(0,time)
   ewise f(i)=%if(i<=%nreg,x(i)*e(time),e(time)**2-sigsq)
   compute v=v+%outerxx(f)
end do time
*
compute v=v*%nobs
compute vinv=inv(v)
*
dim %hstats(%nreg+1)
compute %hjoint=0.0,%hstats=%const(0.0),s=%const(0.0)
do time=startl,endl
   if .not.smpll(time)
      next
   compute x=%eqnxvector(0,time)
   ewise f(i)=%if(i<=%nreg,x(i)*e(time),e(time)**2-sigsq)
   compute s=s+f
   ewise %hstats(i)=%hstats(i)+s(i)**2
   compute %hjoint=%hjoint+%qform(vinv,s)
end do time
ewise %hstats(i)=%hstats(i)/v(i,i)

if .not.print
   return

report(action=define,hlabels=||'Test','Statistic','P-Value'||)
report(atrow=1,atcol=1) 'Joint' %hjoint HTestCrit(%hjoint,%nreg+1)
report(atrow=2,atcol=1) 'Variance' %hstats(%nreg+1) HTestCrit(%hstats(%nreg+1),1)
do i=1,%nreg
   report(atrow=2+i,atcol=1) %eqnreglabels(0)(i) %hstats(i) HTestCrit(%hstats(i),1)
end do i
report(action=format,atcol=2,tocol=2,width=10)
report(action=format,atcol=3,tocol=3,picture='#.##')
report(action=show,window='Hansen Stability Test')
end stabtest

