/* This program is an adaptation of code originally written by R. Carter Hill at LSU. The program deals with Example 4.9.4 in Greene, 4th edition, p. 157. */ options charcode nocenter linesize=80; data expmle; input obs inc edu; cards; 1 20.5 12 2 31.5 16 3 47.7 18 4 26.2 16 5 44.0 12 6 8.28 12 7 30.8 16 8 17.2 12 9 19.9 10 10 9.96 12 11 55.8 16 12 25.2 20 13 29.0 12 14 85.5 16 15 15.1 10 16 28.5 18 17 21.4 16 18 17.7 20 19 6.42 12 20 84.9 16 ; proc iml; start mle; use expmle; read all into x var{edu}; read all into y var{inc}; print , "Maximum Likelihood Estimation of Gamma Distribution", "Greene, 4th ed., p. 157"; /* Starting values for beta and rho */ theta = {-4,3}; crit = 1; t = nrow(x); result = j(10,7,0); /* Use Newton's method to maximize likelihood function. Note the following SAS functions: lgamma(rho)=log of gamma function evaluated at rho; digamma(rho)=first derivative of log of gamma function evaluated at rho; trigamma(rho) = second derivative of log of gamma function evaluated at rho. */ do iter = 1 to 20 while (crit > 1.0e-10); beta = theta[1,1]; rho = theta[2,1]; logl = -rho*log(beta+x) - lgamma(rho) - y/(beta+x) + (rho-1)*log(y); logl = sum(logl); /* Calculate gradient vector */ g1 = sum (-rho/(beta+x) + y/((beta + x)##2)); g2 = sum (-log(beta+x) - digamma(rho) + log(y) ); g = g1//g2; /* Form the Hessian */ h11 = sum ((rho/(beta+x)##2) - 2*y/((beta+x)##3)); h22 = -t * trigamma(rho); h12 = -sum (1/(beta+x)); h = (h11||h12)//(h12`||h22); /* The Newton two step */ db = -inv(h)*g; thetanew = theta + db; crit = sqrt(ssq(thetanew - theta)); /* Place the result of this iteration in the result matrix */ result[iter,] = iter||g1||g2||(theta`)||crit||logL; theta = thetanew; end; cov = -inv(H); * cov from Hessian; i11 = -rho*sum(1/(beta+x)##2); info = - ((i11||h12)//(h12||h22)); covmle = inv(info); * cov from Information matrix; * recall that x is fixed; * and E(y|x)=rho*(beta+x); g1 = (-rho/(beta+x) + y/((beta + x)##2)); g2 = (-log(beta+x) - digamma(rho) + log(y) ); gmat = g1||g2; bhhh = gmat`*gmat; covbh3 = inv(bhhh); * cov from BHHH; cnames = {iter,g1,g2,beta,rho,crit,logl}; print "result",result [colname=cnames]; pnames = {beta,rho}; print , "The Maximum Likelihood Estimates-Unrestricted Model", theta [rowname=pnames]; print , "Hessian", H [rowname=pnames colname=pnames]; print , "Asymptotic Covariance Matrix-From Hessian", cov [rowname=pnames colname=pnames]; print , "Information Matrix", info [rowname=pnames colname=pnames]; print , "Asymptotic Covariance Matrix-From Information Matrix", covmle [rowname=pnames colname=pnames]; print , "BHHH estimate of information matrix", bhhh [rowname=pnames colname=pnames]; print , "Asymptotic Covariance Matrix-From bhhh", covbh3 [rowname=pnames colname=pnames]; Wald = (rho-1)*inv(cov[2,2])*(rho-1); * Wald test; critval = cinv(.95,1); * calculates the 95th percentile of chi-square 1; pval = 1 - probchi(wald,1); * calculates the probability value of Wald; print "Results of Wald test Using Hessian" Wald critval pval; Wald = (rho-1)*inv(covmle[2,2])*(rho-1); * Wald test; critval = cinv(.95,1); * calculates the 95th percentile of chi-square 1; pval = 1 - probchi(wald,1); * calculates the probability value of Wald; print "Results of Wald test Using Info Matrix" Wald critval pval; Wald = (rho-1)*inv(covbh3[2,2])*(rho-1); * Wald test; critval = cinv(.95,1); * calculates the 95th percentile of chi-square 1; pval = 1 - probchi(wald,1); * calculates the probability value of Wald; print "Results of Wald test Using BHHH" Wald critval pval; loglu = logl; *start rmle; print , "Maximum Likelihood Estimation of Exponential Dist", "Greene, p. 157"; beta = 10; crit = 1; t = nrow(x); result = j(10,5,0); do iter = 1 to 20 while (crit > 1.0e-10); logl = -log(beta+x) - y/(beta+x); logl = sum(logl); g = sum (-1/(beta+x) + y/((beta + x)##2)); h = sum ((1/(beta+x)##2) - 2*y/((beta+x)##3)); db = -inv(h)*g; betanew = beta + db; crit = sqrt(ssq(betanew - beta)); result[iter,] = iter||g||beta||crit||logL; beta = betanew; end; cov = -inv(h); cnames = {iter,g,beta,crit,logl}; print ,result [colname=cnames]; pnames = {beta}; print , "The Maximum Likelihood Estimates-Restricted Model", beta [rowname=pnames]; print , "Asymptotic Covariance Matrix-From Hessian", cov [rowname=pnames colname=pnames]; print , "Hessian", H [rowname=pnames colname=pnames]; print "Evaluation of Unrestricted Model Log-likelihood at RMLE",; rho = 1; logl = -rho*log(beta+x) - lgamma(rho) - y/(beta+x) + (rho-1)*log(y); logl = sum(logl); loglr = logl; /* Gradient evaluated at restricted MLE estimates */ g1 = sum (-rho/(beta+x) + y/((beta + x)##2)); g2 = sum (-log(beta+x) - digamma(rho) + log(y) ); g = g1//g2; /* Hessian evaluated at restricted MLE estimates */ h11 = sum ((rho/(beta+x)##2) - 2*y/((beta+x)##3)); h22 = -t * trigamma(rho); h12 = -sum (1/(beta+x)); h = (h11||h12)//(h12`||h22); Print , "restricted Log-likelihood = " logl, "gradient",g,"hessian",h; /*Information Matrix evaluated at restricted MLE estimates */ i11 = -rho*sum(1/(beta+x)##2); info = - ((i11||h12)//(h12||h22)); covrmle = inv(info); LM = g`*(-inv(h))*g; * LM test; critval = cinv(.95,1); pval = 1 - probchi(LM,1); print "Results of LM test Using Hessian" LM critval pval; LM = g`*covrmle[2,2])*g; * LM test; critval = cinv(.95,1); pval = 1 - probchi(LM,1); print "Results of LM test Using Info Matrix" LM critval pval; g1 = (-rho/(beta+x) + y/((beta + x)##2)); g2 = (-log(beta+x) - digamma(rho) + log(y) ); gmat = g1||g2; bhhh = gmat`*gmat; covbh3r = inv(bhhh); LM = g`*covbh3r[2,2]*g; * LM test; critval = cinv(.95,1); pval = 1 - probchi(LM,1); print "Results of LM test Using BHHH" LM critval pval; LR = -2*(loglr-loglu); * Likelihood Ratio test; pval = 1 - probchi(LR,1); print "Results of LR test" LR critval pval; finish; run mle; quit;