7z'>f "f|%GRAD.M: calculates gradient and Hessian % %INPUTS: % b1: starting values for the estimates; % y,X: data and size of data % nobs:size of data % %OUTPUT: % g:gradient of LLF % H: Hessian of LLF function [g, H]=grad(y,X,bsv,nobs) gg=zeros(nobs,2); %--------------------------compute gradient-------------------- for i=1:nobs gg(i,:)=y(i)*X(i,:)-X(i,:)*exp(bsv'*X(i,:)'); end; g=sum(gg)/nobs; %--------------------------compute Hessian-------------------- H=zeros(2,2); for i=1:nobs H=H+(X(i,:)'*X(i,:))*exp(bsv'*X(i,:)'); end; H=-H/nobs; %LOGL.M: calculates log-likelihood %for the Poisson regression % %INPUTS: % b1: starting values for the estimates; % y,X: data and size of data % nobs:size of data % %OUTPUT: % LL: value of log-likelihood function % g:gradient of LLF % H: Hessian of LLF function [LL, g, H] = logl(b1,X,y,nobs,varargin); %--------------------------compute likelihood-------------------- ll=zeros(nobs,1); for i=1:nobs ll(i)=y(i)*b1'*X(i,:)'-exp(b1'*X(i,:)'); end; LL=-sum(ll); %negative sign is used to get maximization, not minimization %--------------------------compute gradient-------------------- if nargout>1 for i=1:nobs gg(i,:)=y(i)*X(i,:)-X(i,:)*exp(b1'*X(i,:)'); end; g=-sum(gg)/nobs; %negative sign is used to get maximization, end; %--------------------------compute Hessian-------------------- if nargout>1 H=zeros(2,2); for i=1:nobs H=H+(X(i,:)'*X(i,:))*exp(b1'*X(i,:)'); end; H=H/nobs; %positive sign is used to get maximization, end; %----------------------------------------------------------- % Part A: Newton-Raphson algorithm (manual implementation) %----------------------------------------------------------- clear all; clc; fid=fopen('out.out','w'); fprintf(fid,'*******************************************************************************\n'); fprintf(fid,' "out.out" generated by script "newton.m" \n'); fprintf(fid,'*******************************************************************************\n'); fprintf(fid,'main routine "newton.m" requires two functions:\n'); fprintf(fid,'- grad.m: computes gradient and hessian for the LLF in the manual NR algorithm in part "A"\n'); fprintf(fid,'- logl.m: computes LLF, gradient and hessian (all with the opposite signs for the maximization) in part "B"\n'); status=fclose(fid); %---------- 1. Generating data ------------------------ nobs=200; X=[ones(nobs,1) unifrnd(0,1,nobs,1)]; b_true=[1.5;-0.5]; y=zeros(nobs,1); for i=1:nobs lambda=exp(b_true'*X(i,:)'); y=poissrnd(lambda,nobs,1); end %---------- 2. Specifying estimation options ------------------------ sv=[0.001;0.001]; %starting values opt.iter=500; %maximum number of iteration allowed opt.conver=1.0000e-8; %convergence criteria opt.count=0; %actual iterations counter %---------- 3. Manual N-R procedure ------------------------ b1=sv; for i=1:opt.iter [g, H]=grad(y,X,b1,nobs); %compute gradient and Hessian for LLF bnew1=b1+inv(-H)*g'; %compute the new value of parameter estimates tol=abs((bnew1-b1)./b1); if all(tol|t| | [95 perc. Conf.Interval]\n'); fprintf(fid,'------------|---------|------------|-----------|----------|--------------------------\n'); for i=1:k fprintf(fid,'%s', char(names(i))); fprintf(fid,' | %5.5f | %5.5f | %5.5f | %5.5f | [%5.5f %5.5f] \n',b1(i),se_b(i),tstat(i),pval(i),CI.a(i),CI.b(i)); end fprintf(fid,'------------------------------------------------------------------------------------\n\n\n'); status=fclose(fid); %---------------------------------------------------------------------------------------------- %---------- 6. Part B. Estimating regression via optimization toolbox ------------------------ %---------------------------------------------------------------------------------------------- b2=sv; options=optimset('fminunc'); options=optimset('LargeScale', 'on' , ... 'GradObj' , 'on' , ... 'Hessian' , 'on' , ... 'Display','off' , ... 'MaxIter' , 500 , ... 'TolX', 1.0000e-8 , ... 'TolFun' , 1.0000e-8, ... 'DerivativeCheck' , 'off' , ... 'Diagnostics' , 'off' , ... 'MaxFunEvals', 2000); [bnew2,fval,exitflag,output,grad,hessian]=fminunc('logl',b2,options,X,y,nobs); %---------- 7. Regression statistics ------------------------ k=size(X,2); % get number of regressors == columns in X se_b=sqrt(diag(inv(hessian)/nobs)); % s.e.(b) tstat=bnew2./se_b; % t-statistics pval=2*(1-tcdf(abs(tstat),nobs-k)); % p-value CI.a=bnew2-se_b; % confidence interval CI.b=bnew2+se_b; %---------- 8. Printing output ------------------------ fid=fopen('out.out','a'); names={'const ' 'X2 '}; fprintf(fid,'\n'); fprintf(fid,'*******************************************************************************\n'); fprintf(fid,' Part B: Optimization toolbox results \n'); fprintf(fid,' using "fminunc" function \n'); fprintf(fid,'*******************************************************************************\n'); fprintf(fid,'Convergence achieved in = %6.0f iterations\n',output.iterations); fprintf(fid,'SSE = %5.5f\n',SSE); fprintf(fid,'-------------------------------------------------------------------------------------\n'); fprintf(fid,'Variable | Coef. | Std.Err. | t | P>|t| | [95 perc. Conf.Interval]\n'); fprintf(fid,'------------|---------|------------|-----------|----------|--------------------------\n'); for i=1:k fprintf(fid,'%s', char(names(i))); fprintf(fid,' | %5.5f | %5.5f | %5.5f | %5.5f | [%5.5f %5.5f] \n',bnew2(i),se_b(i),tstat(i),pval(i),CI.a(i),CI.b(i)); end fprintf(fid,'------------------------------------------------------------------------------------\n\n\n'); fprintf(fid,'---------------END OF FILE ---------------------\n'); status=fclose(fid); b1 bnew23O]{U;*ɞfI]W%u3;hB!mX_Tq¨ h4;".ceۊbd݁}0r\ڧ"/<~ r #] 9z