/*
** Lesson 12.1: GMM Estimation of a Gamma Distribution
** See Greene [1999], Example 4.26 and Example 11.4
*/
use gpe2;
output file=gpe\output12.1 reset;

load data[21,2]=gpe\yed20.txt;
x=data[2:21,1]/10;      @ income: data scaling may help @ 

call reset;
_nlopt=0;               @ it is a minimization problem @
_method=5;
_iter=100;

_b={3,1};             @ initial values of parameters @
_hacv=1;                @ hetero consistent covariance @
                        @ assuming serially uncorrelated @
call estimate(&gmmqw,x);
                        @ using the results of previous estimator @
_b=__b;                 @ for initial value of parameters and @
gmmw=invpd(gmmv(x,_b)); @ for computing the weight matrix @
call estimate(&gmmqw,x); 
call gmmout(x,__b);     @ print GMM output @
/*
_b=__b;
call estimate(&gmmq,x);
call gmmout(x,__b);     @ print GMM output @
*/
end;

/*
User-defined moments equations, must be named mf
based on gamma distribution: b[1]=rho, b[2]=lambda
*/
proc mf(x,b);
    local n,m;
    n=rows(x);
    m=zeros(n,4);
    m[.,1]=x-b[1]/b[2];
    m[.,2]=x^2-b[1]*(b[1]+1)/(b[2]^2);
    m[.,3]=ln(x)-gradp(&lngamma,b[1])+ln(b[2]);
    m[.,4]=1/x-b[2]/(b[1]-1);
    retp(m);
endp;

/*
Log of gamma distribution function
*/
fn lngamma(x)=ln(gamma(x));

#include gpe\gmm.gpe;