options memory=30 crt limwarn=0 ;
?
? ========================= SIMCMLEHS2 ===============================
?
?  One Variable Time Series Process Using Panel Data.
?  Use simulated data to investigate testing methodology.
?  Estimator: CMLE with homoskedastic errors.
?  DGPs have heteroskedastic errors.
?
?  Bronwyn H. Hall - Hall and Mairesse 2001 (for Munich).
?
supres @covoc @smpl ;

const nob 200 t 12 ndraw 10 ;
set nobt = nob*t ;
set t1 = t-1 ;
smpl 1 nobt ;
list dylist dy79-dy89 ;
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 ;
				 
do imod = 1 to 1 ;
  date startime ;
  do isim = 1 to ndraw ;
  smpl 1 nobt ;
	input simhsprcs ;    ? Simulate the 8 models.
	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(corr,silent) dylist ;
      unmake @mean mean79-mean89 ;
	  unmake @var var79-var89 ;
	dot 79-89 ;
      dy. = dy.-mean. ;
    enddot ;
	
?   Compute likelihood for rho=1.
    set sigsq = (var79+var80+var81+var82+var83+var84
	            +var85+var86+var87+var88+var89)/t1 ;
	set sigsq = sigsq*(nob-1)/nob ;   ? Use ML Estimate of sigsq.
    set logl = -(nob*t1/2)*(log(2*3.14159)+1+log(sigsq)) ;
	set logl0(isim) = logl ;

?   Estimate rho.	
	set srho = 1 ;
	set rho = 0.5 ;
	cmlecx dylist rho srho var svar logl conv ;
?	print rho srho var svar logl conv ;
   	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 ; set crit=@miss ;
    set reject(isim) = crit<=.05 ;
	set lrreject(isim) = crit2<=.05 ;
  enddo ;
  smpl 1 ndraw ;
?  print rhoest srhoest lltest lrtest logl0 ;
  dot rhoest srhoest vareps svare lltest reject logl0 lrtest lrreject
      ifconv ;  
    unmake . . ;
  enddot ;
  msd (pair) rhoest srhoest lltest reject vareps svare logl0 lrtest
             lrreject ifconv ;
  dot rhoest srhoest vareps svare lltest reject logl0 lrtest lrreject
      ifconv ;  
    . = @miss ;
  enddot ;
  
  date endtime ;
  set time = (endtime-startime)/60 ;
  print (format="( ' Elapsed time in minutes:',f6.1)") time ;
enddo ;