options memory=40 crt limwarn=0 ;
?
? ========================= SIMCMLEHSHS ===============================
?
?  One Variable Time Series Process Using Panel Data.
?  Use simulated data to investigate testing methodology.
?  Estimator: CMLE with heteroskedastic errors.
?  DGPs have heteroskedastic errors.
?  Bronwyn H. Hall - Hall and Mairesse 2001 (for Munich).
?
supres @covoc @smpl ;
const nob 200 t 12 ndraw 1 ;
smpl 1 ndraw ;
list stats ymean yvar rhoest srhoest vareps svare lltest reject
     logl0 lrtest lrreject ;
dot stats ;	   		 ? Initialize simulation statistics.
  . = @miss ;
enddot ;
  
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 ;

date startime ;
do imod = 1 to 8 ;
  do isim = 1 to ndraw ;
  smpl 1 nobt ;
  input simhsprcs ;   ? Simulate the 8 models; stationary and heteroskedastic errors.
  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(corr,silent) dylist ;
      unmake @mean mean79-mean89 ;
    dot 79-89 ;
      dy. = dy.-mean. ;
    enddot ;
	
?   Compute likelihood for rho=1.
    set srho = -1 ;
	set rho = 1 ;
	cmlechsx dylist rho srho var svar logl conv ;
?	print rho srho var svar logl conv ;
	set logl0(isim) = logl ;
	smpl 1 nob ;

?   Maximize likelihood with respect to rho and sigsq. 
	set srho = 1 ;
	param rho .5 ;
	cmlechsx dylist rho srho var svar logl conv ;
?	print rho srho var svar logl 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 ;
  dot stats ;	   
    unmake . . ;
  enddot ;
  msd (pair) stats ;
  dot stats ;
    . = @miss ;
  enddot ;

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