gmm.m (2591B)
1 function [beta,se,Omega]=gmm(beta0,X,uname,maxit,param,tolerance,verbose); 2 3 % function [beta,se,Omega]=gmm(beta0,X,uname,maxit,param,tolerance,verbose); 4 % 5 % GMM, following Ogaki notes (paper M056). 12/2/97 6 % The idea is very straightforward. We have a collection of q moment 7 % restrictions: E z'e = 0. Define u=f(X,beta)=z'e (can be nonlinear). 8 % 9 % GMM involves choosing beta to minimize the quadratic form 10 % 11 % J(beta) = ubar'*W*ubar where ubar is the sample mean of the u 12 % 13 % and W is some "weighting matrix". 14 % 15 % Define Omega = Eu'u. Then the optimal weighting matrix is inv(Omega). 16 % In practice, we calculate Omega by iterating: Assume W=I, get beta, get 17 % a sample Omega, set W=inv(omega), and iterate until convergence. 18 % 19 % Then, se(beta) = sqrt(diag(1/T inv(Gamma'inv(Omega)*Gamma))), 20 % where Gamma is the gradient of u wrt beta. 21 % 22 % Then Hansen J test of overidentifying restrictions is T*J ~ X^2(q-p). 23 % 24 % 25 % beta0 = initial guess 26 % X = collection of all data needed to compute moments 27 % uname = function that takes beta0,X,param and returns moments u 28 % maxit = max iterations 29 % param = parameters to be passed to u function. 30 % tolerance= minimum percentage gain in Jfct before stopping 31 % verbose= 0 if quiet 32 % 33 % Note well: uname must return a qxT matrix of u's, not qx1! 34 % (This is required for the computation of Omega). 35 % -- In OLS, se's correspond exactly to White Robust se's. 36 37 38 options=foptions; 39 if exist('verbose')==1; options(1)=verbose; end; 40 if exist('tolerance')~=1; tolerance=.001; end; 41 if exist('param')~=1; param=1; end; 42 if exist('maxit')==1; options(14)=maxit; end; 43 44 options(3)=1e-6; % default = 1e-4 45 beta=beta0; 46 getu=['u=' uname '(beta,X,param);']; 47 eval(getu); 48 [q T]=size(u); 49 if q>T; disp 'error q>T'; abc; end; 50 p=size(beta,1); 51 52 % First, use W=I 53 W=eye(q); 54 gain=100; 55 Jold=1e+10; 56 57 % Now iterate until tolerance achieved 58 while gain>tolerance; 59 beta=fmins('Jfct',beta,options,[],X,W,uname,param) 60 eval(getu); 61 Omega=1/T*u*u'; 62 Jnew=Jfct(beta,X,W,uname,param); 63 gain=abs((Jnew-Jold)/Jold); 64 Jold=Jnew; 65 Wold=W; 66 W=inv(Omega); 67 end; 68 W=Wold; 69 70 Gt=zeros(q,p); 71 for t=1:T; 72 Gt=Gt+gradp(uname,beta,[],X(t,:),param); 73 end; 74 G=1/T*Gt; 75 76 vcv=1/(T-p)*inv(G'*inv(Omega)*G); 77 se=sqrt(diag(vcv)); 78 tstat=beta./se; 79 80 JStat=T*Jnew; 81 dof=q-size(beta,1); 82 pval=1-chi2cdf(JStat,dof); 83 84 fprintf('J: %12.4f\n',Jnew); 85 fprintf('JStat: %5.2f Pval: %5.3f DoF: %5.3f\n',[JStat pval dof]); 86 disp ' Beta S.E. t-stat '; 87 disp '--------------------------------------'; 88 cshow(' ',[beta se tstat],'%12.5f');