7z¼¯'ªît %"Ȉ>_function h = halton(lgth,prim); % This function constructs a Halton sequence of length lgth based on a % user-provided prime number. The function returns a lgthx1 vector % containing the Halton sequence. The values are 'well-spaced' % realizations from the (0,1) interval. They can be used directly in % simulation or transformed into other random variable realizations. %create first part of the sequence k = floor(log(lgth+1)/log(prim)); phi = 0; i = 1; while i <= k; x = phi; j = 1; while j < prim; y = phi + (j/prim^i); x = [x; y]; j = j + 1; end; phi = x; i = i + 1; end; % create remainder of sequence x = phi; j = 1; while (j < prim) & (length(x) < lgth + 1); y = phi + (j/prim^i); x = [x; y]; j = j + 1; end; % return sequence of length lgth h = x(2:lgth+1); % Maximum Likelihood estimation of the Mixed Logit Model %----------------------------------------------------------- clear all; clc; diary out.out; warning off; %----------------------------------------------------------- % 1. Single Data Set %----------------------------------------------------------- N=40; % Number of observations J=4; % Number of categories k=4; % number of parameters to estimate NNC=2; % number of random coefficients R=20; % number of draws for the Halton sequence var_names={' b1 ' ' b3 ' ' b2 ' ' s2 ' ' b4 ' ' s4 '}; X=unifrnd(0,3,N,J*k); % generetate X b1=ones(N,1); % generate b() b2=3+0.75*randn(N,1); % b3=ones(N,1); % b4=3+0.25*randn(N,1); % b=[b1 b2 b3 b4]; % e=log(-log(unifrnd(0,1,N,J))); % generate disturbance % compute utility U = zeros(N,J); for i=1:k U = U + b(i)* X(:,(i-1)*J+1:i*J)+e; end; % compute Y [C,I] = max(U,[],2); % get indices of the alternatives with the highest utility Y=zeros(N,J); for i=1:N Y(i,I(i))=1; end % compute halton sequences for each individual errors=zeros(R,N*NNC); for i=1:N*NNC h = halton(R,3)'; errors(:,i)=norminv(h,0,1); end; % starting values b=[0.25;0.25;0.25;0.25]; % starting values for parameters q=[0.25;0.25]; % starting values for pdf(NC) b=[b;q]; FC=[1;3]; % positions of of fixed coefficents NC=[2;4]; % positions of of normal coefficents % set-up optimization options options=optimset('fminunc'); options=optimset('LargeScale', 'on' , ... 'GradObj' , 'off' , ... 'Hessian' , 'off' , ... 'Display','off' , ... 'MaxIter' , 100 , ... 'TolX', 1.0000e-6 , ... 'TolFun' , 1.0000e-6, ... 'DerivativeCheck' , 'off' , ... 'Diagnostics' , 'off' , ... 'MaxFunEvals', 2000); % evaluate the model [b_mxl,fval,exitflag,output,grad,hessian]=fminunc('mxl_like',b,options,Y,X,FC,NC,errors); b_mxl(4)=abs(b_mxl(4)); %adjust for possible non-positive estimate b_mxl(6)=abs(b_mxl(6)); %adjust for possible non-positive estimate % estimation statistics stder=diag(sqrt(inv(hessian))); tstat=b_mxl./stder; pval=2*(1-tcdf(abs(tstat),N-k-NNC)); % printing output clc; fprintf ('**********************************************************************\n'); fprintf (' 1. Mixed Logit Estimation. Single Data Set\n'); fprintf ('**********************************************************************\n\n'); fprintf ('Number of Hamilton draws: R = %6.0f\n',R); fprintf ('Number of observations: N = %6.0f\n\n',N); fprintf('-----------------Coefficients on explanatory variables----------------\n'); fprintf(' Coeff. | St.err. | t-stat. | p-val.\n'); fprintf('----------------------------------------------------------------------\n'); for i=1:size(b_mxl) fprintf('%s', char(var_names(i))); fprintf (' %6.6f %6.6f %6.6f %6.6f\n',b_mxl(i),stder(i),tstat(i),pval(i)); end; fprintf('----------------------------------------------------------------------\n'); %----------------------------------------------------------- % 2. Multiple Data Sets %----------------------------------------------------------- N=40; % Number of observations NNC=2; % number of random coefficients R=100; % number of draws for the Halton sequence M=2; % number of replications b_mxl1=zeros(M,k+2); %specify matrix with estimated coefficients over the replications X=unifrnd(0,3,N,J*k); % generetate X % compute halton sequences for each individual errors=zeros(R,N*NNC); for i=1:N*NNC h = halton(R,3)'; errors(:,i)=norminv(h,0,1); end; for loops=1:M loops b1=ones(N,1); % generate b() b2=3+0.75*randn(N,1); b3=ones(N,1); b4=3+0.25*randn(N,1); b=[b1 b2 b3 b4]; e=log(-log(unifrnd(0,1,N,J))); % generate disturbance % compute utility U = zeros(N,J); for i=1:k U = U + b(i)* X(:,(i-1)*J+1:i*J)+e; end; % compute Y [C,I] = max(U,[],2); Y=zeros(N,J); for i=1:N Y(i,I(i))=1; end % starting values b=[0.25;0.25;0.25;0.25;0.25;0.25]; % starting values for parameters FC=[1;3]; % positions of of fixed coefficents NC=[2;4]; % positions of of normal coefficents % set-up optimization options options=optimset('fminunc'); options=optimset('LargeScale', 'on' , ... 'GradObj' , 'off' , ... 'Hessian' , 'off' , ... 'Display','off' , ... 'MaxIter' , 100 , ... 'TolX', 1.0000e-6 , ... 'TolFun' , 1.0000e-6, ... 'DerivativeCheck' , 'off' , ... 'Diagnostics' , 'off' , ... 'MaxFunEvals', 2000); % evaluate the model b_mxl1(loops,:)=fminunc('mxl_like',b,options,Y,X,FC,NC,errors)'; b_mxl1(4)=abs(b_mxl1(4)); %adjust for possible non-positive estimate b_mxl1(6)=abs(b_mxl1(6)); %adjust for possible non-positive estimate end b_mxl_loop=mean(b_mxl1); % estimation statistics stder_loop=std(b_mxl1); tstat_loop=b_mxl_loop./stder_loop; pval_loop=2*(1-tcdf(abs(tstat_loop),N-k-NNC)); % printing output fprintf ('**********************************************************************\n'); fprintf (' 2. Mixed Logit Estimation. Multiple Data Sets\n'); fprintf ('**********************************************************************\n\n'); fprintf ('Number of Hamilton draws: R = %6.0f\n',R); fprintf ('Number of observations: N = %6.0f\n',N); fprintf ('Number of simulations: M = %6.0f\n\n',M); fprintf('-----------------Coefficients on explanatory variables----------------\n'); fprintf(' Coeff. | St.err. | t-stat. | p-val.\n'); fprintf('----------------------------------------------------------------------\n'); for i=1:size(b_mxl) fprintf('%s', char(var_names(i))); fprintf (' %6.6f %6.6f %6.6f %6.6f\n',b_mxl_loop(i),stder_loop(i),tstat_loop(i),pval_loop(i)); end; fprintf('----------------------------------------------------------------------\n'); diary out; figure(1) for i=1:size(b_mxl1,2) subplot(3,2,i); hist(b_mxl1(:,i)); title(['beta',num2str(i)]); end; warning on;function ll = mxl_like(b,Y,X,FC,NC,errors); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % A function to simulate the likelihood function for a specific version of % the mixed logit model that contains fixed and independent random % coefficients. % % Inputs include: % b: sx1 coefficient vector % Y: nxj matrix of observed outcomes % X: nx(j*k) matrix of explanatory variables % FC: NFCx1 vector indicating positions of fixed coefficents % NC: NNCx1 vector indicating positions of normal coefficients % errors: Rx(n*NNC) matrix of standard normal realizations % % The function returns the scalar simulated (negative) log-likelihood %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % establish dimensions for model [R junk] = size(errors); [n j] = size(Y); [n jk] = size(X); k = jk/j; s = length(b); NFC = length(FC); NNC = length(NC); % initialize probability vector p0 = zeros(n,j); % Add in fixed coefficients V = zeros(n,j); for h = 1:NFC; g = FC(h); V = V + b(h)* X(:,(g-1)*j+1:g*j); h = h + 1; end; % Add in normal coefficients and calculate simulated probabilities % by looping over r=1,...R. for r = 1:R; err = errors(r,:); err = reshape(err,n,NNC); EV = V; for h = 1:NNC; g = NC(h); mu = ones(n,j)*b(NFC+(2*h)-1); dev = kron(ones(1,j),err(:,h))*b(NFC+(2*h)); EV = EV + (mu + dev).*X(:,(g-1)*j+1:g*j); h = h + 1; end; EV = exp(EV); den = kron(sum(EV')',ones(1,j)); p0 = p0 + EV./den; r = r + 1; end; lnp = log(p0/R); ll = -sum(sum((Y.*lnp)')); 3®O¥$Í]ñԡΊÁ7jN¸Vüœœã†ïq ¯E4lñë:§Xò‚ÔN/aNQ¡: ãŒæJóém”µÄEøJS©XçžbìB÷W\CKN¨£M½k¡ îòò@Ojª½7;âìàT&×Å(sÏ*P ù§€¥" ~ #] €¤ 6…D