' THIS PROGRAM PERFORMS 2 LM SPECIFICATION TESTS FOR THE LOGIT MODEL:
' (1) An LM test of the null that the null hypothesis that the errors of a Logit model are homoskedastic, against the alternative that they are heteroskedastic, being a function of the "r" variables in Z. We will call this the LM-HET test
'NOTE: Z MUST NOT INCLUDE AN INTERCEPT - See Davidson & MacKinnon (2004, p.464)
' (2) An LM test of the null that the model is Logistic, against a more flexible functional form. We will call this the LM_RESET test.
' The Illlustrative code below involves a model with with 2 covariates and an intercept. In this example the "Z" variables that give rise to the potential heteroskedasticity are taken to be the 2 covariates in the model itself . This can be amended to suit - see the comments in the code below.
'The notation and construction of the test statistics follows the exposition in Davidson & MacKinnon (2004, pp.464-465).
' References:
' 1. Davidson, R. & J. G. MacKinnon (1984), "Convenient specification tests for logit and probit models", Journal of Econometrics, 25, 241-262.
' 2. Davidson, R. & J. G. MacKinnon (2004), "Econometric Theory and Methods", Oxford University Press, Oxford
'=============================================================
' This version of theode written by David Giles, Dept. of Economics, University of Victoria, October 2011
'=============================================================
' Adjust the sample size to suit
!n=5000
smpl 1 !n
series bf
series sf
series v
series lhs
series one=1
scalar r
' Modify the covariates in the group "xx" and in the Lrobit model, as appropriate.
' Also modify the hetero. variates in the group "zz" as appropriate, and set "r" equal to the number of variables in "zz". (NOTE: "zz" must NOT include an intercept variable)
' Note that Y is the binary (zero-one) dependent variable
group xx one x1 x2
matrix x=@convert(xx)
r=2
group zz x1 x2
matrix z=@convert(zz)
' Note that in the next line, the Logit model is estimated with het-consistent standard errors - see the "h" option. This could be changed if you want.
equation eq1.logit(h) y xx
fit(i) xbeta
bf=1/(one+exp(-xbeta))
sf=exp(xbeta)/(one+exp(xbeta))^2
v=bf*(one-bf)
lhs=(y-bf)/@sqrt(v)
series scalefacx=sf/@sqrt(v)
series scalefacz= -xbeta*scalefacx
' The line immediately above was wrong in the earlier version of this code
stom(scalefacx, scalexv)
stom(scalefacz, scalezv)
matrix r1m=@scale(x, scalexv)
matrix r2m=@scale(z, scalezv)
mtos(r1m,r1)
mtos(r2m,r2)
group bigr r1 r2
' Fit the auxiliary regression for the homoskedasticity test:
equation eq2.ls lhs bigr
fit fitted
series f2=fitted*fitted
scalar lm_het=@sum(f2)
' The homoskedasticity test statistic is asymptotically Chi Square with "r" degrees of freedom, under the null of homoskedasticity, where "r" is the number of variables involved in the (alternative) heteroskedasticity. So calculate the (asymptotic) p-value:
scalar het_p_val=1-@cchisq(lm_het,r)
' Fit the auxiliary regression for the functional form (RESET) test:
series r3=scalefacx*(xbeta)^2
group rhs r3 r1
' The coeffcient to be tested is that of the first regerssor - see the ordering of the above group
equation eq2.ls lhs rhs
scalar lm_reset=@tstats(1)
' The functional form test statistic is asymptotically standard normal. So calculate the (asymptotic) p-value:
scalar reset_p_val=2*(1-@cnorm(abs(lm_reset)))