options limwarn=0 memory=20 crt ;
supres smpl ;
?
? ========================= SIMIPS0 ===============================
?
?  One Variable Time Series Process - Simulation Results. 
?  All processes; homoskedastic errors
?  Im, Pesaran, & Shim test for unit roots.
?  p=0 only.
?  This version uses simulation of stationary processes.
?  BCHM paper Avril 1998.
?  Revised September 2001.
?

const nob 2000 t 12 ndraw 100 ;
set nobt = nob*t ;
smpl 1 nobt ;
random (seedin=29345) e ; ? Initialize seed.
?
?  Set up id and trend.
?
trend(pstart=1,period=t) trend ;
select trend=1 ;
  trend id ;
select trend>1 ;
  id = id(-1) ;
select 1 ;
  dummy trend yrdum ;
  
date startime ;
do imod = 1 to 8 ;
?
?   Initialize result storage
?
  dot not t ;
    mform(nrow=ndraw,ncol=1) ips. ;
    mform(nrow=ndraw,ncol=1) pval. ;
  enddot ;
  mform(nrow=ndraw,ncol=1) at ;
  mform(nrow=ndraw,ncol=1) dr ;

  smpl 1 nob ;
  dot f0 f0t ;
    ad. = 0 ;
  enddot ;
  adfid = 0 ;
  alpha = 0 ;
  delta = 0 ;

?     Loop over ndraw simulations.

  do isim = 1 to ndraw ;
    smpl 1 nobt ;
    input simprocs ;    ? Simulate the 8 stationary models.
    set j = 0 ;
?    
?      Remove year-specific means.
?
    olsq (silent) y yrdum ;
    copy @res y ;
?
?      Loop over individual unit root tests (t time periods per obs).
? 
    do ind = 1 to nobt by t ;        
      set ibeg = ind ;
      set iend = ind+t-1 ;
      smpl ibeg iend ;             
        set j = j+1 ;
        unit (notrend,nows,silent) y ;? Test with no trend.
        set adf0(j) = @tabdf(2,1) ;
        unit (nows,silent) y ;        ? Unit root test with trend.
        set adf0t(j) = @tabdf(2,1) ;
        set adfid(j) = id(ind) ;
        set alpha(j) = @tabdf(4,4) ;
?       if alpha(j)>1000 ; then ; set outlier = ind ;
        set delta(j) = @tabdf(6,4) ;
      enddo ;

?     Now compute the IPS tests by averaging the individual tests.

      smpl 1 nob ;
      msd (silent) adf0;        

      set zp0not = sqrt(j)*(@mean(1)+1.506)/sqrt(1.04) ;?  p=0 or 2 or 3 t=12
?     title 'Tests with augmenting lags = 0' ;
      cdf(noprint,lowtail,normal) zp0not crit ;
      set pvalnot(isim) = crit<=.05 ;
      set ipsnot(isim) = zp0not ;

?     Do the same thing for the test with a trend. 

      msd (silent) alpha delta ;
      set at(isim) = @mean(1) ;
      set dt(isim) = @mean(2) ;
      msd (silent) adf0t ; 
      set zp0t = sqrt(j)*(@mean(1)+2.166)/sqrt(1.08) ;

?     title 'Test with augmenting lags 0, plus time trend' ;
      cdf(noprint,lowtail,normal) zp0t crit ;
      set pvalt(isim) = crit<=.05 ;
      set ipst(isim) = zp0t ;
    enddo ;

  smpl 1 ndraw ;

  dot not t ;
    unmake ips. zp0. ;
    unmake pval. pval0. ;
    msd zp0. pval0. ;
  enddot ;
  unmake at alpha ;
  unmake dt delta ;
  msd alpha delta ;
enddo ;

date endtime ;
set time = (endtime-startime)/60 ;
print (format="( ' Elapsed time in minutes:',f7.2)") time ;