*
* @RegSTRTest( options )
* does an LM test for linearity vs an alternative of smooth transition
* based upon lags of a known threshold variable in the preceding linear
* regression. This should be executed immediately after a LINREG.
*
* Options:
*   THRESHOLD=series for threshold variable [dependent variable of regression]
*   D        =delay on threshold variable [1]
*   TITLE    =title for output ["Test for STR in regression"]
*   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 lagged threshold
* variable.
*
*   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 LSTR or ESTR should cause rejection of linearity and
* rejection of H12. H12 is the appropriate statistic if ESTR is the main
* hypothesis of interest. If H12 is rejected, but H03 is accepted, that
* points towards an ESTR model, as opposed to an LSTR.
*
* References:
*   The general idea (of testing interactions between regressors and terms
*   from a Taylor expansion of the transition variable) is from
*
*     Luukkonen, Saikkonen and Terasvirta(1988), "Testing Linearity against
*     smooth transition autoregressive models", Biometrika vol 75, 491-499.
*
*   The labeling and interpretation of the hypotheses is from
*
*     Terasvirta(1994), "Specification, Estimation and Evaluation of Smooth
*     Transition Autoregressive Models", JASA, vol 89, no. 425, pp 208-218.
*
* Revision Schedule:
*  04/2013 Written by Tom Doan, Estima
*
procedure RegSTRTest
*
option series  threshold
option integer d      1
option switch  print  1
option string  title
option series  weights
*
* Y is the local copy of the threshold variable
*
local series   y
local rect[series] trans
*
local integer 	i j k actualdf test
local vect     rss fstat signif
local rect     s stest
local string   ltitle
local report   sreport
*
* Generate the cross products of the regressors with the first three
* powers of the threshold variable, so these will involve 2nd to 4th
* powers of the threshold variable.
*
if %defined(threshold)
   set y = threshold
else
   set y = %eqndepvar(0){0}
*
compute k=%nreg
dim trans(k,3)
do j=1,3
   do i=1,k
      set trans(i,j) = %eqnxvector(0,t)(i)*(y{d}^j)
   end do i
end do j
*
* Compute the cross product matrix of the regressors, the transformed
* variables and the residuals.
*
cmom(weights=weights)
# %rlfromtable(%eqntable(0)) trans %resids
dim rss(4) fstat(5) signif(5)
*
* Sweep out standard regressors
*
compute s=%sweeptop(%cmom,k)
compute rss(1)=s(%ncmom,%ncmom)
*
* Test for linear dependence among the regressors and the first power
* interactions. We will have this if the transition variable is one of
* the regressors or a linear combination of regressors). Figure out how
* many degrees of freedom we actually have in each batch of auxiliary
* regressors.
*
compute actualdf=k
compute stest=%decomp(%xsubmat(%cmom,1,2*k,1,2*k))
do test=1,k
   if stest(k+test,k+test)==0.0
      compute actualdf=actualdf-1
end do test
*
* 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*k+1,(j+1)*k))
   compute rss(j+1)=s(%ncmom,%ncmom)
end do j
*
* Compute the f-stats for the sequence of hypotheses.
*
compute %ndf=%nobs-k
do j=1,3
   compute fstat(j+1) =((rss(j)-rss(j+1))/actualdf)/(rss(j+1)/(%ndf-j*actualdf))
   compute signif(j+1)=%ftest(fstat(j+1),k,%ndf-j*actualdf)
end do j
compute fstat(5)=(rss(1)-rss(3))/(2*actualdf)/(rss(3)/(%ndf-2*actualdf))
compute signif(5)=%ftest(fstat(5),2*actualdf,%ndf-2*actualdf)
*
* And the overall test.
*
compute %cdstat=(rss(1)-rss(4))/(3*actualdf)/(rss(4)/(%ndf-3*actualdf))
compute %signif=%ftest(%cdstat,3*actualdf,%ndf-3*actualdf)

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Test for STR in regression"

report(use=sreport,action=define,title=ltitle)
report(use=sreport,atrow=1,atcol=1,span) ltitle
report(use=sreport,atrow=2,actol=2)  "Threshold" %l(threshold)
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
