7z'3b!d%LOGL.M: calculates log-likelihood for the multimomial logit model % %INPUTS: % b1: starting values for the estimates; % y,X: data % %OUTPUT: % LL: value of log-likelihood function % g:gradient of LLF function [LL] = logl(b1,X,Y,varargin); %--------------------------compute likelihood-------------------- [nobs,j]=size(Y); k=size(b1,1); V=zeros(nobs,j); for i=1:k V=V+b1(i)*X(:,(i-1)*j+1:i*j); end; den=sum(exp(V'))'; den1=kron(den,ones(1,j)); P=exp(V)./den1; LL = -sum(sum(Y.*log(P))); g1=X'*(Y - P); for i=1:k g(i,:)=sum(sum(g1((i-1)*j+1:i*j,:))); end; g=-g/nobs; % Maximum Likelihood estimation of the Multinomial Logit model % Main script "mlogit.m" also requires function "logl.m" %----------------------------------------------------------- %---------- 1. Specify output file ------------------------ clear all; clc; fid=fopen('output.out','w'); fprintf(fid,'*******************************************************************************\n'); fprintf(fid,' "output.out" generated by script "mlogit.m" \n'); fprintf(fid,'*******************************************************************************\n\n'); fprintf(fid,'main routine "mlogit.m" requires the following function:\n'); fprintf(fid,'- logl.m: computes LLF and gradient (both with the opposite signs for the maximization)\n\n'); status=fclose(fid); %---------- 2. Generating data ------------------------ nobs=200; J=4; X1=2*rand(nobs,J); X2=rand(nobs,J); X=[X1 X2]; b=[1;2.5]; EV=-log(-log(rand(nobs,4))); U=X1*b(1)+X2*b(2)+EV; [c,index]=max(U'); index=index'; Y=zeros(nobs,J); zz=1; for i=1:nobs Y(i,index(i,1))=1; end; %-------------3. Likelihood estimation ----------------- k=size(b,1); sv=[1.5;1.5]; options=optimset('fminunc'); options=optimset('LargeScale', 'on' , ... 'GradObj' , 'off' , ... 'Hessian' , 'off' , ... 'Display','off' , ... 'MaxIter' , 500 , ... 'TolX', 1.0000e-8 , ... 'TolFun' , 1.0000e-8, ... 'DerivativeCheck' , 'off' , ... 'Diagnostics' , 'off' , ... 'MaxFunEvals', 2000); [bnew,fval,exitflag,output,grad,hessian]=fminunc('logl',sv,options,X,Y); %-------------4. Estimation statistics ----------------- % 4.a Get it from Hessian ----------------- se_b=sqrt(diag(inv(hessian)/nobs)); tstat=bnew./se_b; % t-statistics pval=2*(1-tcdf(abs(tstat),nobs-k)); % p-values CI.a=bnew-se_b; % confidence intervals CI.b=bnew+se_b; % 4.a Get it score*score' ----------------- se_b1=sqrt(diag(grad*grad')); tstat1=bnew./se_b1; % t-statistics pval1=2*(1-tcdf(abs(tstat1),nobs-k)); % p-values CI.a1=bnew-se_b1; % confidence intervals CI.b1=bnew+se_b1; %---------- 5. printing estimation output ------------------------- fid=fopen('output.out','a'); names={'X1 ' 'X2 '}; fprintf(fid,'*******************************************************************************\n'); fprintf(fid,' Maximum Likelihood estimation of the Multinomial Logit model \n'); fprintf(fid,'*******************************************************************************\n'); fprintf(fid,'Number of observations = %6.0f\n',nobs); fprintf(fid,'Number of parameters to estimate = %5.0f\n',k); fprintf(fid,'True parameters: b1= %2.5f b2= %2.5f\n',b(1),b(2)); 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',bnew(i),se_b(i),tstat(i),pval(i),CI.a(i),CI.b(i)); end fprintf(fid,'------------------------------------------------------------------------------------\n\n'); fprintf(fid,'In the table above all the statistics are calculated using numerical Hessian\n'); fprintf(fid,'When the product of scores is used to get t-statistics, the results are as follows:\n'); 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',bnew(i),se_b1(i),tstat1(i),pval1(i),CI.a1(i),CI.b1(i)); end fprintf(fid,'------------------------------------------------------------------------------------\n\n'); status=fclose(fid); %------------- End of file ----------------- 31O= i1pEm;?9xh!GI!DwaJ\s4x ? opewgzBr 9jz ] #] o ^