new;
 output file = test.out reset ;
library optmum;

N = 12;

beta0 = .5; beta1=3 ; sigma1 = 1; beta2=0 ; sigma2 = 4; rho=.8 ;

x=seqa(1,1,N)/N.*rndn(N,1) ; /* x is the "regressor," could be anything, but if not centered around 0,
the estimates will not be well identified. As for OLS you would like ``variation in the regressor"
which here mean that we want to observe a good spread between 0s and 1s */

z=ones(N,1)+2*x+5*rndn(N,1) ; /* z is another regressor for linear
relation ---not orthogonal to x*/

Xmat=ones(N,1)~x~z ;

Nsim =1 ; "Number of simulations:" Nsim ; /* create a matrix to
hold the estimated values */

results = zeros(Nsim,3) ;


/*starting values*/

let sv = 1 1 1 1;

let svR = 1 1 1 ;

estim=0 ;
do while estim < Nsim ;
estim=estim+1 ;

u = sigma1*rndn(N,1);

epsilon=sigma2*rndn(N,1) ;


/*generate a linear relation */

y=beta0+beta1*x+beta2*z+epsilon ;

__output=0;

{b_ml,lmax,g,retcode} = optmum(&logl,sv);

 {b_mlR,lmaxR,g,retcode}=optmum(&loglR,svR) ;

__output=1 ;

/*results[estim,.]=b_ml' ;*/


endo ; " " ;

/*results of last estimation*/

 "calculate Hessian"; h=hessp(&logl,b_ml);

 /*this is minus the hessian because we minimize -L*/

 "variance matrix"; vmat=inv(h);
 "OLS variance matrix" ; OLSvmat=inv(xmat'xmat).*b_ml[4] ;
 /*vmat; */

 "" ;

"beta2-hat, t-value" b_ml[3] b_ml[3]/sqrt(OLSvmat[3,3])

" " ; "P-values:" ;

pval1=2*cdftc(abs(b_ml[3])/sqrt(vmat[3,3]),N-4) ;

 "t-test:" ;  pval1;

/* pval2=cdfchic(PUT LR TEST HERE,1) ;

"LR test" ; pval2;

*/

 pval3=cdfchic((b_ml[3])^2/vmat[3,3],1) ;

"Wald test" ; pval3 ;

/*restricted ML-estimator with restriction*/

b_LM=zeros(4,1) ;

b_LM[1]=b_mlR[1] ; b_LM[2]=b_mlR[2] ; b_LM[4]=b_mlR[3] ;

LMgrad=gradp(&logL,b_LM);

h_LM=hessp(&logL,b_LM);

var_LM=inv(h_LM) ;

LMtest=LMgrad*var_LM*LMgrad' ;

pval4=cdfchic(LMtest,1) ;

 "LM test" ; pval4 ;


"***";
" " ;

proc logl(b);
  local r,L,muvec,sigma2,beta0,beta1,beta2,beta3,beta4, i;

beta1=b[1] ;
beta2=b[2] ;
beta3=b[3] ;
sigma2=b[4] ;

L=0 ;
i=0 ;
do while i<N ;
  i=i+1 ;
    L = L  -0.5*ln(sigma2^2) - 0.5*(y[i]-beta1-beta2*x[i]-beta3*z[i])^2/sigma2^2;
endo ;
  retp(-L); /*optmum minimizes*/
endp;

proc loglR(b);
  local r,L,muvec,sigma2,beta0,beta1,beta2,beta3,beta4, i;

beta1=b[1] ; beta2=b[2] ; sigma2=b[3] ;

L=0 ; i=0 ; do while i<N ;
  i=i+1 ;
    L = L  -0.5*ln(sigma2^2) - 0.5*(y[i]-beta1-beta2*x[i])^2/sigma2^2;

endo ;
  retp(-L); /*optmum minimizes*/
endp; output off ;