nlls5.m (2136B)
1 function [u,beta,se,tstat]=nlls5(beta0,X,mname,options,param,meth,uname); 2 3 % nlls5.m 4/17/96 (from phinlls.m in ~/rad). 4 % 5 % nlls5 is just nlls.m, updated for Version 5 of matlab (which uses 6 % optimset and fminsearch/fminunc) as new optimization routines. 7 % 8 % A generic NLLS routine, to be adapted to whichever problem is at hand 9 % 10 % beta0=initial guess 11 % X = [y x] 12 % 13 % Model: y=m(beta,x)+u 14 % 15 % The function uufct calls m(beta,x) to compute u'*u. 16 % mname is a user-defined function (e.g. 'mfct') giving 17 % m(beta,x). 18 % 19 % options=passed (default is foptions) 20 % options=foptions; 21 % options(1)=verbose; 22 % options(14)=maxit; 23 % options(2)=tol; 24 % options(3)=tol; 25 26 % meth= "u" for fminu or "s" for fmins (default is fmins); 27 % uname = optional u'*u function; default is ~/procs/uufct (standard) 28 29 if exist('options')~=1; options=optimset('fminsearch'); end; 30 if isempty(options); options=optimset('fminsearch'); end; 31 if exist('param')~=1; param=[]; end; 32 if exist('meth')~=1; meth='s'; end; 33 if exist('uname')~=1; uname='uufct'; end; 34 35 if optimget(options,'MaxIter')>0; 36 if meth=='s'; 37 beta=fminsearch(uname,beta0,options,X); 38 else; 39 beta=fminunc(uname,beta0,options,X); 40 end; 41 else; 42 disp 'Just evaluating R2 at beta0...'; 43 beta=beta0; % if just evaluating R2 44 end; 45 46 % Now calculate standard errors using the Gradients 47 [N col]=size(X); 48 y=X(:,1); 49 x=X(:,2:col); 50 k=length(beta); 51 if isempty(param); 52 M=gradp(mname,beta,[],x); 53 eval(['u=y-' mname '(beta,x);']); 54 else; 55 M=gradp(mname,beta,[],x,param); 56 eval(['u=y-' mname '(beta,x,param);']); 57 end; 58 sigma2=u'*u/(N-k); 59 mminv=inv(M'*M); 60 vcv=sigma2*mminv; 61 se=sqrt(diag(vcv)); 62 tstat=beta./se; 63 64 % White robust vcv 65 f=mult(M,u); 66 Sw=f'*f/(N-k); 67 robvcv=N*mminv*Sw*mminv; 68 robse=sqrt(diag(robvcv)); 69 robt =beta./robse; 70 71 ybar=1/N*(ones(1,N)*y); 72 rsq=1-(u'*u)/(y'*y - N*ybar^2); 73 74 fprintf('R-Squared: %7.4f\n',rsq); 75 fprintf('uu*100: %7.4f\n',u'*u*100); 76 disp ' Beta S.E. t-stat WhiteSE WhiteT'; 77 disp '------------------------------------------------------------'; 78 cshow(' ',[beta se tstat robse robt],'%12.6f'); 79 80