nllsnw.m (1757B)
1 function [u,beta,se,tstat]=nllsnw(beta0,X,mname,verbose,maxit); 2 3 % nllsnw.m 4/17/96 (from phinllsnw.m in ~/rad). 4 % A generic NLLSNW routine, to be adapted to whichever problem is at hand 5 % 6 % beta0=initial guess 7 % X = [y x] 8 % 9 % Model: y=m(beta,x)+u 10 % 11 % The function uufct calls m(beta,x) to compute u'*u. 12 % mname is a user-defined function (e.g. 'mfct') giving 13 % m(beta,x). 14 % 15 % 2/25/97 -- the 'nw' extension -- the robust errors are newey-west serial 16 % correlation robust. 17 18 options=foptions; 19 if exist('verbose')==1; options(1)=verbose; end; 20 if exist('maxit')==1; options(14)=maxit; end; 21 22 beta=fmins('uufct',beta0,options,[],X,mname); 23 24 % Now calculate standard errors using the Gradients 25 [N col]=size(X); 26 y=X(:,1); 27 x=X(:,2:col); 28 k=length(beta); 29 M=gradp(mname,beta,[],x); 30 eval(['u=y-' mname '(beta,x);']); 31 sigma2=u'*u/(N-k); 32 mminv=inv(M'*M); 33 vcv=sigma2*mminv; 34 se=sqrt(diag(vcv)); 35 tstat=beta./se; 36 37 % White robust vcv 38 %f=mult(M,u); 39 %Sw=f'*f/(N-k); 40 %robvcv=N*mminv*Sw*mminv; 41 %robse=sqrt(diag(robvcv)); 42 %robt =beta./robse; 43 44 % Newey-West robust vcv (from lsnw.m) 45 G=2*floor(N^(1/4))+1; 46 B=0; 47 for j=0:G; 48 lam=0; 49 for t=j+1:N; 50 st= (M(t,:))'*u(t); 51 stj=(M(t-j,:))'*u(t-j); 52 lam = lam +st*stj'; 53 end % t 54 lam = lam/N; 55 B=B+(1-j/(G+1))*(lam+lam'); 56 end %j 57 B=(N/(N-k))*B; 58 robvcv =mminv*(N*B)*mminv; 59 robse = sqrt(diag(robvcv)); 60 robt = beta./robse; 61 62 ybar=1/N*(ones(1,N)*y); 63 rsq=1-(u'*u)/(y'*y - N*ybar^2); 64 65 fprintf('R-Squared: %7.4f\n',rsq); 66 fprintf('uu*100: %7.4f\n',u'*u*100); 67 disp ' Beta S.E. t-stat NewWstSE NewWstT'; 68 disp '------------------------------------------------------------'; 69 cshow(' ',[beta se tstat robse robt],'%12.6f'); 70 71