options memory=40 crt double limwarn=0 ;
?
? ========================= SIMCMLE ===============================
?
?  One Variable Time Series Process Using Panel Data.
?  Use simulated data to investigate testing methodology.
?  CMLE with homoskedastic disturbances (Kruiniger 1999).
?  DGP can be homoskedastic (simprocs) or heteroskedastic (simhsprcs).
?  Bronwyn H. Hall - Hall and Mairesse 2001 (for Munich).
?
supres @covoc @smpl ;
const nob 200 t 12 ndraw 1000 ;
list stats ymean yvar rhoest srhoest lltest reject vareps svare
           logl0 lrtest lrreject ifconv ;
smpl 1 ndraw ;          ?  Initialize the simulation statistics.
dot stats ;
  . = @miss ;
enddot ;

set nobt = nob*t ;
set t1 = t-1 ;
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 trend>0 ;
  dummy trend yrdum ;
list dylist dy79-dy89 ;	   ? List of differenced series.

do imod =1 to 8 ;          ? Simulate all 8 models (see labels below).
  date startime ;
  do isim = 1 to ndraw ;
    smpl 1 nobt ;

?   Choose between simprocs and simhsprcs below, depending on which process
?   you wish to simulate.
    input simprocs ;     ? Simulate the 8 models using stationary homoskedastic DGPs. 
?    input simhsprcs ;   ? Simulate the 8 models using stationary heteroskedastic DGPs. 

	msd (silent) y ;
	set ymean(isim) = @mean(1) ;
	set yvar(isim) = @var(1) ;
	select trend>1 ;
      dy = y-y(-1) ;
      mmake dymat dy ;
      mform (nrow=t1,ncol=nob) dymat ;
      mat dymat=dymat' ;
    smpl 1 nob ;
      unmake dymat dylist ;
      msd(silent) dylist ;
      unmake @mean mean79-mean89 ;
	  unmake @var var79-var89 ;
    dot 79-89 ;
      dy. = dy.-mean. ;
    enddot ;

?   Compute likelihood for rho=1 using ML Estimate of sigsq,
?   which is easily computed from directly from the data.

?   set srho = -1 ;  ? Alternative method uses MLE with rho==1.
?   set rho = 1 ;	 ? This method is slower, but gives same result.
?   cmlecx dylist rho srho var svar logl conv ; 

	set sigsq = (var79+var80+var81+var82+var83+var84
	            +var85+var86+var87+var88+var89)/t1 ;
	set sigsq = sigsq*(nob-1)/nob ;   
    set logl = -(nob*t1/2)*(log(2*3.14159)+1+log(sigsq)) ;
	set logl0(isim) = logl ;

	set srho = 1 ;
	param rho 0.5 var 1 ;
	set var = sigsq ;
    cmlecx dylist rho srho var svar logl conv ; ? Estimate by CMLE.
	if srho=1 ; then ; srho=@miss ;
	set ifconv(isim) = conv ;
	set vareps(isim) = var ;
	set svare(isim) = svar ;
	set rhoest(isim) = rho ;
	set srhoest(isim) = srho ;
    set test = (rho-1)/srho ;
	set test2 = (logl-logl0(isim))*2 ;
	if abs(rho-1)<.000001 ; then ; do ;
      set test = 0 ; set test2 = 0 ;
	enddo ;
    set lltest(isim) = test ;
	set lrtest(isim) = test2 ;
	if ~miss(test) ; then ; do ;
	  cdf (normal,lowtail) test crit ;
	  cdf (chisq,df=1) test2 crit2 ; 
	  enddo ;
	  else ; do ;
	    set crit = @miss ;
		set crit2 = @miss ;
	  enddo ;
    set reject(isim) = crit<=.05 ;
	set lrreject(isim) = crit2<=.05 ;
  enddo ;
  
  smpl 1 ndraw ;
  dot stats ;
    unmake . . ;
  enddot ;
  msd (pair) stats ;
  dot stats ;
	. = @miss ;
  enddot ;

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