*
*
*
* Estimation of threshold
*
* Control parameters
*
* Fraction of threshold values excluded at each end
*
compute pi=.15
*
* Number of VAR lags
*
compute maxlag=4
*
* Threshold delay
*
compute d = 1
*
* Moving average terms in smoothing threshold variable
*
compute malength = 2
*
*cal(q) 2000
*allocate 2014:1
*


open data FISCAL.xlsx
calendar(q) 2000:03
data(format=xlsx,org=columns) 2000:03 2014:01 gdp spen inf r inves ex tax

* These are the desired start and end of the sample to be used in the
* estimation. However, the transformations use the entire range.

*
compute sstart = 2000:3
compute send   = 2014:1

*******************************************
          system(model=varmodel)
          variables gdp spen inf r
          lags 1 to maxlag
          det constant
          end(system)
********************************************
compute rstart=sstart+maxlag,rend=send
estimate(print,resids=u) rstart rend
*
compute loglr=%logl
compute nvar =%nvar
*
filter(type=lagging,span=malength) gdp / gdpthr
*
* Get a sorted copy of the threshold variable
*
set copy rstart rend = gdpthr{d}
order copy rstart rend
*
compute piskip=fix(pi*%nobs)+(maxlag*nvar+1)
compute pistart=rstart+piskip,piend=rend-piskip
*
compute bestlogl=loglr
*
set lrstats pistart piend = 0.0
do pientry=pistart,piend
   compute thresh=copy(pientry)
   *
   * This allows for the covariance matrix to differ between regimes
   *
   sweep(var=hetero,group=gdpthr{d}<thresh) rstart rend
   # %modeldepvars(varmodel)
   # %rlfromeqn(%modeleqn(varmodel,1))
   *
   * If this is best so far, save it
   *
   if %logl>bestlogl
       compute bestlogl=%logl,bestthresh=thresh
   compute lrstats(pientry)=2.0*(%logl-loglr)
end do pientry
*
disp "Best Threshold Value" bestthresh
compute suplr=2.0*(bestlogl-loglr)
*
* The degrees of freedom of the test count both the number of VAR
* parameters and the number of elements in the second estimated
* covariance matrix.
*
disp "Degrees of freedom" %nregsystem/2+%nvar*(%nvar+1)/2
*
sstats(mean) pistart piend lrstats>>avglr log(lrstats)>>avgloglr
compute explr=exp(avgloglr)
*
* Fixed regressor bootstrap
*
dec vect[series] fiddled(%nvar)
*
compute nboot=500
set avglrboot 1 nboot = 0.0
set explrboot 1 nboot = 0.0
set suplrboot 1 nboot = 0.0
*
do reps=1,nboot
   set lrstats pistart piend = 0.0
   set fiddlef = %ran(1.0)
   do i=1,%nvar
      set fiddled(i) = fiddlef*u(i)
   end do i
   *
   sweep rstart rend
   # fiddled
   # %rlfromeqn(%modeleqn(varmodel,1))
   compute loglr=%logl
   compute bestlogl=loglr
   *
   do pientry=pistart,piend
      compute thresh=copy(pientry)
      *
      * This allows for the covariance matrix to differ between regimes
      *
      sweep(var=hetero,group=gdpthr{d}<thresh) rstart rend
      # fiddled
      # %rlfromeqn(%modeleqn(varmodel,1))
      *
      * If this is best so far, save it
      *
      if %logl>bestlogl
          compute bestlogl=%logl,bestthresh=thresh
      compute lrstats(pientry)=2.0*(%logl-loglr)
   end do pientry
   sstats(mean) pistart piend lrstats>>avglrboot(reps) log(lrstats)>>avgloglr
   compute explrboot(reps)=exp(avgloglr)
   compute suplrboot(reps)=2.0*(bestlogl-loglr)
end do reps
sstats(mean) 1 nboot (explrboot>explr)>>exppvalue $
                     (suplrboot>suplr)>>suppvalue $
                     (avglrboot>avglr)>>avgpvalue
*
report(action=define)
report(atrow=1,atcol=1) "" "Statistic" "P-value"
report(atrow=2,atcol=1) "Sup" suplr suppvalue
report(atrow=3,atcol=1) "Avg" avglr avgpvalue
report(atrow=4,atcol=1) "Exp" explr exppvalue
report(atcol=2,tocol=2,action=format,width=8)
report(atcol=3,tocol=3,action=format,picture="*.###")
report(action=show)