尤度比検定、Wald検定、ＬＭ検定　(プログラム)

// LR test, Wald test, LM test after Maximum likelihood estimation
#include<oxstd.h>
#include<oxfloat.h> // to use M_2PI
#import<maximize>   // to use maximization package

static decl s_mX, s_vY; // static variable to use in the functions
// Unconstrained Loglikelihood Function
fLoglik(const vP, const adFunc, const avScore, const amHess)
{
decl ck,cnobs;decl dsum,dsig2;decl vb,ve;
 cnobs=rows(s_vY);ck=rows(vP);
 vb=vP[0:ck-2];dsig2=vP[ck-1];ve=s_vY-s_mX*vb;
 dsum=-sumsqrc(ve)/(2*dsig2);
 adFunc[0]=-0.5*cnobs*log(M_2PI*dsig2)+dsum;
if(avScore) // provides analytical 1st derivative
 {
 (avScore[0])[0:ck-2]=s_mX'*ve/dsig2;
 (avScore[0])[ck-1]=-0.5*cnobs/dsig2+sumsqrc(ve)/(2*dsig2^2);
 }
if(amHess) // provides expected value of hessian matrix
 {
 (amHess[0])=zeros(ck,ck);
 (amHess[0])[0:ck-2][0:ck-2]=-s_mX'*s_mX/dsig2;
 (amHess[0])[ck-1][ck-1]=-cnobs/(2*dsig2^2);
 }
 return 1;
}
// Constrained Loglikelihood Function
fLoglikC(const vP, const adFunc, const avScore, const amHess)
{
decl ck,cnobs;decl dsum,dsig2;decl vb,ve;
 cnobs=rows(s_vY);ck=rows(vP);
 vb=vP[0:ck-2]|1-vP[2]-vP[3];dsig2=vP[ck-1];
 ve=s_vY-s_mX*vb;dsum=-sumsqrc(ve)/(2*dsig2);
 adFunc[0]=-0.5*cnobs*log(M_2PI*dsig2)+dsum;
return 1;
}
main()
{
decl ck,cnobs; // scalar (count, number of ...)
decl dfunc,dlm,ds2,dssr,dloglr,dloglu,dlrt; // scalar (double)
decl mhess,mi,msigma,mx; // matrix
decl vp,vpvalue,vse,vscore,vstat,vzvalue; // vector
decl file;
 // Read file
 file = fopen("nerlove.asc");fscan(file,"%#m",145,5,&mx);
 fclose(file); mx=log(mx);
//
 println("Maximum Likelihood Esimation using BFGS method.");
 s_mX=mx[][1:2]~mx[][4]~mx[][3];  // X matrix: independent variable
 s_mX=ones(145,1)~s_mX; //          + constant term
 s_vY=mx[][0]; // Y vector: dependent variable
 cnobs=rows(s_mX); // number of observations
 ck=columns(s_mX)+1; // number of parameters (including constant term & sigma^2)
 vp=<0;0;0;0;0;10>; // initial values

println("*************************************************");
println("Unconstrained Model");
println("*************************************************");
//
// Unconstrained maximum likelihood estimation
// y=b1+b2*x2+b3*x3+b4*x4+b5*x5
//
 MaxBFGS(fLoglik,&vp,&dfunc,0,TRUE);
 dloglu=dfunc; // log L (unconstrained)
 vscore=zeros(ck,1); // score vector (d logL/d theta)
 fLoglik(vp,&dfunc,&vscore,&mhess);
 mi=-mhess;msigma=invert(mi);
 vse=sqrt(diagonal(msigma))';
 vzvalue=vp./vse; // z-value for H0: parameter=0
 vpvalue=tailn(fabs(vzvalue)); 
 println("%r",{"beta1","beta2","beta3","beta4","beta5","sigma^2"},
         "%c",{"est. coeff.", "std. err.","z-value", "p-value"},
         "%10.6",vp~vse~vzvalue~vpvalue);
 // Wald Test
 decl dg,vg,dwald; // g & G
  dg=vp[2]+vp[3]+vp[4]-1;vg=<0,0,1,1,1,0>;
  dwald=dg*invert(vg*msigma*vg')*dg;

println("*************************************************");
println("Constrained Model (B3+B4+B5=1)");
println("*************************************************");//
// Constrained maximum likelihood estimation
// y=b1+b2*x2+b3*x3+b4*x4+(1-b3-b4)*x5
//
 decl vpc;
 vpc=vp[0:3]|vp[5];
 MaxBFGS(fLoglikC,&vpc,&dfunc,0,TRUE);
 dloglr=dfunc; // log L (constrained)
 Num2Derivative(fLoglikC,vpc,&mhess); // Numerical Hessian Matrix
 mi=-mhess;msigma=invert(mi); 
 vse=sqrt(diagonal(msigma))';
 vzvalue=vpc./vse; // z-value for H0: parameter=0
 vpvalue=tailn(fabs(vzvalue)); 
 println("%r",{"beta1","beta2","beta3","beta4","sigma^2"},
         "%c",{"est. coeff.", "std. err.","z-value", "p-value"},
         "%10.6",vpc~vse~vzvalue~vpvalue);
 dlrt=-2*(dloglr-dloglu);
 // LM Test
  vp=vpc[0:3]|1.0-vpc[2]-vpc[3]| vpc[4];
  fLoglik(vp,&dfunc,&vscore,&mhess);
  mi=-mhess;msigma=invert(mi);
  dlm=vscore'*msigma*vscore;
  vstat= dlm | dlrt | dwald;vpvalue=tailchi(vstat,1);
println("*************************************************");
println("TEST H0:B3+B4+B5=1 vs H1:B3+B4+B5 =! 1");
println("*************************************************");

 println("%r",{"LM test", "Likelihood ratio test", "Wald test"},
         "%c",{"Statistic", "p-value"},
         "%10.6",vstat~vpvalue);
}