*
* @StarTest( options ) y start end
* does an LM test for linearity vs an alternative of LSTAR or ESTAR
* (smooth transition autoregression) in a univariate autoregression.
*
* Options:
*   D     = delay on threshold variable [1]
*   P     = lags in the autoregression [1]
*   TITLE = title for output ["Test for STAR in series ..."]
*   WEIGHTS=series of observation weights [equally weighted]
*   [PRINT]/NOPRINT
*
* Variables Defined:
*   %CDSTAT  overall linearity LM statistic
*   %SIGNIF	 significance level of %CDSTAT
*
* Test Statistics:
*
*   Linearity is the joint test of all the interaction terms between the
*   regressors and the powers (1, 2 and 3) of the lags.
*
*   H01 is a test of the first order interaction terms only
*   H02 is a test of the second order interaction terms only
*   H03 is a test of the third order interaction terms only
*   H12 is a test of the first and second order interactions terms only
*
* Either an LSTAR or ESTAR should cause rejection of linearity and
* rejection of H12. H12 is the appropriate statistic if ESTAR is the
* main hypothesis of interest. If, H12 is rejected, but H03 is accepted,
* that points towards an ESTAR model, as opposed to an LSTAR.
*
* References:
*  Terasvirta(1994), "Specification, Estimation and Evaluation
*  of Smooth Transition Autoregressive Models", JASA, vol 89, no. 425, pp
*  208-218.
*
* Revision Schedule:
*  07/2005 Written by Tom Doan, Estima
*  03/2008 Updated to add comments
*  08/2010 Use local report. Reformat slightly. Add TITLE and PRINT options.
*  11/2010 Add output for H12 testing the 1st and 2nd order terms.
*  11/2011 Correct the labeling on the hypotheses
*
procedure StarTest y start end
type series 	y
type integer 	start end
*
option integer d      1
option integer p      1
option switch  print  1
option string  title
option series  weights
*
local rect[series] trans

local integer 	i j
local vect 		rss fstat signif
local rect 		s
local string 	ltitle
local report   sreport

if .not.%defined(y) {
   display "Syntax: @STARTest(options) >>y<< start end"
   return
}
*
* Generate the cross products of the p lags of y with the first three
* powers of the threshold variable, so these will involve 2nd to 4th
* powers of y's.
*
dim trans(p,3)
do j=1,3
   do i=1,p
      set trans(i,j) = y{i}*(y{d}^j)
   end do i
end do j
*
* Compute the cross product matrix of the standard AR(p) regressors and
* the transformed variables.
*
cmom(weights=weights) start end
# constant y{1 to p} trans y
dim rss(4) fstat(5) signif(5)
*
* Sweep out standard regressors
*
compute s=%sweeplist(%cmom,%seq(1,p+1))
compute rss(1)=s(%ncmom,%ncmom)
*
* Sweep out the lags, starting from the second order to the fourth. Save
* the sum of squared residuals after each.
*
do j=1,3
   compute s=%sweeplist(s,%seq(j*p+2,(j+1)*p+1))
   compute rss(j+1)=s(%ncmom,%ncmom)
end do j
*
* Compute the f-stats for the sequence of hypotheses.
*
compute %ndf=%nobs-p-1
do j=1,3
   compute fstat(j+1) =((rss(j)-rss(j+1))/p)/(rss(j+1)/(%ndf-j*p))
   compute signif(j+1)=%ftest(fstat(j+1),p,%ndf-j*p)
end do j
compute fstat(5)=(rss(1)-rss(3))/(2*p)/(rss(3)/(%ndf-2*p))
compute signif(5)=%ftest(fstat(5),2*p,%ndf-2*p)
*
* And the overall test.
*
compute %cdstat=(rss(1)-rss(4))/(3*p)/(rss(4)/(%ndf-3*p))
compute %signif=%ftest(%cdstat,3*p,%ndf-3*p)

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Test for STAR in series "+%l(y)

report(use=sreport,action=define,title=ltitle)
report(use=sreport,atrow=1,atcol=1,span) ltitle
report(use=sreport,atrow=2,atcol=1)  "AR length" p
report(use=sreport,atrow=3,atcol=1)  "Delay"     d
report(use=sreport,atrow=5,atcol=1)  "Test" "F-stat" "Signif"
report(use=sreport,atrow=6,atcol=1)  "Linearity"  %cdstat %signif
report(use=sreport,atrow=7,atcol=1)  "H01" fstat(2) signif(2)
report(use=sreport,atrow=8,atcol=1)  "H02" fstat(3) signif(3)
report(use=sreport,atrow=9,atcol=1)  "H03" fstat(4) signif(4)
report(use=sreport,atrow=10,atcol=1) "H12" fstat(5) signif(5)
report(use=sreport,action=format,atcol=2,width=9)
report(use=sreport,action=format,atcol=3,width=6)
if print
   report(use=sreport,action=show)
end



