nlls.m (1928B)
1 function [u,beta,se,tstat]=nlls(beta0,X,mname,options,param,meth,uname); 2 3 % nlls.m 4/17/96 (from phinlls.m in ~/rad). 4 % A generic NLLS 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 % options=passed (default is foptions) 16 % options=foptions; 17 % options(1)=verbose; 18 % options(14)=maxit; 19 % options(2)=tol; 20 % options(3)=tol; 21 22 % meth= "u" for fminu or "s" for fmins (default is fmins); 23 % uname = optional u'*u function; default is ~/procs/uufct (standard) 24 25 if exist('options')~=1; options=foptions; end; 26 if isempty(options); options=foptions; end; 27 if exist('param')~=1; param=[]; end; 28 if exist('meth')~=1; meth='s'; end; 29 if exist('uname')~=1; uname='uufct'; end; 30 31 if options(14)>0; 32 if meth=='s'; 33 beta=fmins(uname,beta0,options,[],X,mname,param); 34 else; 35 beta=fminu(uname,beta0,options,[],X,mname,param); 36 end; 37 else; 38 beta=beta0; % if just evaluating R2 39 end; 40 41 % Now calculate standard errors using the Gradients 42 [N col]=size(X); 43 y=X(:,1); 44 x=X(:,2:col); 45 k=length(beta); 46 if isempty(param); 47 M=gradp(mname,beta,[],x); 48 eval(['u=y-' mname '(beta,x);']); 49 else; 50 M=gradp(mname,beta,[],x,param); 51 eval(['u=y-' mname '(beta,x,param);']); 52 end; 53 sigma2=u'*u/(N-k); 54 mminv=inv(M'*M); 55 vcv=sigma2*mminv; 56 se=sqrt(diag(vcv)); 57 tstat=beta./se; 58 59 % White robust vcv 60 f=mult(M,u); 61 Sw=f'*f/(N-k); 62 robvcv=N*mminv*Sw*mminv; 63 robse=sqrt(diag(robvcv)); 64 robt =beta./robse; 65 66 ybar=1/N*(ones(1,N)*y); 67 rsq=1-(u'*u)/(y'*y - N*ybar^2); 68 69 fprintf('R-Squared: %7.4f\n',rsq); 70 fprintf('uu*100: %7.4f\n',u'*u*100); 71 disp ' Beta S.E. t-stat WhiteSE WhiteT'; 72 disp '------------------------------------------------------------'; 73 cshow(' ',[beta se tstat robse robt],'%12.6f'); 74 75