'uses of gradients
'hetero-robust Gauss-Newton regression
'see Davidson-MacKinnon 1993, chapter 11.6
'version 4
'last revised 10/27/2000 h.

'change path to program path
%path = @runpath
cd "{%path}"

'create workfile
workfile gallant1 u 1 30

'read data from plain text file
read(b2) amstat.dat t y x1 x2 x3

'starting values
c(1) = -0.04866
c(2) =  1.03884
c(3) = -0.73792
c(4) = -0.51362

'estimate restricted model 
'replicate Gallant (1987) fig 5a (p.35)
equation eq1.ls(c=1e-9,m=1000,showopts,deriv=aa) y = c(1)*x1 + c(2)*x2 + c(4)*exp(c(3)*x3)
show eq1.output

'get residuals
eq1.makeresid res

'get derivative series
eq1.makederivs(n=grd1)

'step 1: regress test variable on derivs and get resids
equation eq2.ls @trend grd1
eq2.makeresids zres

'step 2: weight step 1 resids by restricted resids
series zu = zres*res

'step 3: run robust Gauss-Newton regression (11.69)
eq2.ls 1 zu

'step 4: compute test statistic n*R^2_u
scalar lmstat2 = eq2.@regobs - eq2.@ssr

'show results in table
table out
out(1,1) = "Omitted variable test based on heteroskedasticity robust Gauss-Newton regression:"
out(2,1) = "n*R^2_u = "
out(2,2) = lmstat2
out(2,3) = "p-value = "
out(2,4) = 1 - @cchisq(lmstat2,1)

'one-variable case can be computed by "regular" Gauss-Newton regression with White standard errors
equation eq2.ls(h) res @trend grd1
show eq2.output

'robust t-statistic in test regression
!tstat = eq2.@tstats(1)
'degrees of freedom
!df = eq2.@regobs - eq2.@ncoef
out(3,1) = "robust t-stat = "
out(3,2) = !tstat
out(3,3) = "p-value = "
out(3,4) = 2*(1-@ctdist(@abs(!tstat),!df))
show out

