	% This program computes the random coefficeints discrete choice model
%described
% the "A Research Assistant's Guide to Discrete Choice Models of
%Demand," NBER technical
% paper #221, and "Measuring Market Power in the Ready-to-Eat Cereal
% Industry," NBER WP
% #6387.

	% Written by Aviv Nevo, May 1998.  See http://www.faculty.econ.northwestern.edu/faculty/nevo/supplements/rc_dc_code.htm

	%Eric Rasmusen, Dec. 31, 2005
 %I have added a lot of commands to display matrix size
%I have updated some of the obsolete minimization functions
%I added a diary, a log, with the date you run the program.
 % I wrote this in Matlab 6.
 % I added the separator lines


	%  I also use Bronwyn Hall's, April 2005  modifications for Matlab 7
  %compatibility. See 
%http://emlab.berkeley.edu/users/bhhall/e220c/rc_dc_code.htm

%The inputs to this code are the data matrices ps2.mat and iv.mat

%The output from this program will be results.txt and  mydiary.txt.  . 

disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                    ')


delete mydiary.txt %we don't want any old mydiary.txt file to interfere.
diary mydiary.txt

date = datestr(clock)

	%The "date" command prints out the date, useful for keeping track of the diary.

global invA ns x1 x2 s_jt IV vfull dfull theta1 theti thetj cdid cdindex

	%global ns nmkt nbrn n_inst

	% load data. see description in readme.txt

load ps2.mat
load iv.mat

	 %Use the two files ps2.mat and iv.mat
% http://www.faculty.econ.northwestern.edu/faculty/nevo/supplements/readme.html
% ps2.mat contains the matrices v, demogr, x1, x2, s_jt, id_demo



disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                    ')



ns = 20;       % number of simulated "indviduals" per market %
nmkt = 94;     % number of markets = (# of cities)*(# of quarters)  %
nbrn = 24;     % number of brands per market. if the numebr differs by
%market this requires some "accounting" vector %
 % Default is ns=20, nmkt=94, nbrn=24. Thus we have 24*94=2256 observations
%Set this smaller if you want to have more a dataset you can eyeball better

n_inst = 20   %Number of instruments for price


disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                    ')



IV = [iv(:,2:(n_inst+1)) x1(:,2:(nbrn+1))];

	% The previous line creates a matrix IV from the instruments and the x's.
%20 columns of instruments for price by default, and nbrn brand dummies


IV= IV(1:nmkt*nbrn, :);
s_jt= s_jt(1:nmkt*nbrn, :);
x1=x1(1:nmkt*nbrn, :);
x2= x2(1:nmkt*nbrn, :);

% The previous line  reduces the data matrix, e.g.,  x1, to its first  nmkt*nbrn observations;

disp(['The dimensions of object IV are ' num2str(size(IV)) ])





clear iv


disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                    ')



	% The  vector  below relates each observation to the market it is in %

cdid = kron([1:nmkt]',ones(nbrn,1));

disp(['Object cdid dimensions: ' num2str(size(cdid)) ])


	% the vector below provides for each index the of the last observation   in the data  used here all brands appear in all markets. if this   is not the case the two  vectors, cdid and cdindex, have to be     created in a different fashion but  the rest of the program works fine. 

cdindex = [nbrn:nbrn:nbrn*nmkt]';

disp(['The dimensions of object cdindex are      ' num2str(size(cdindex)) ])

disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	% starting values. zero elements in the following matrix correspond to
% coeff that will not be max over,i.e are fixed at zero.  The defaults are in %comments here

%theta2w=    [0.3302   5.4819         0    0.2037         0;
%             2.4526  15.8935    -1.2000        0    2.6342;
%             0.0163  -0.2506         0    0.0511         0;
%             0.2441   1.2650         0   -0.8091         0];


 theta2w=    [0.3302   5.4819         0    0.2037         0;
              2.4526  15.8935    -1.2000        0    2.6342;
             0.0163  -0.2506         0    0.0511         0;
             0.2441   1.2650         0   -0.8091         0];

	%There are 13 parameters to be estimated, and 7 set to zero by assumption
%These are 4 characteristic parameters and 9 interaction  parameters.
 %Some interactions are set to zero (e.g. Sugar and Incomsquared)


disp('                                                                    ')
disp(['The dimensions of object theta2w are      ' num2str(size(theta2w)) ])
 disp('                                                                    ')

	% create a vector of the non-zero elements in the above matrix, and the
% corresponding row and column indices. this facilitates passing values %
% to the functions below. %

[theti, thetj, theta2]=find(theta2w);

horz=['    mean       sigma      income   income^2    age    child'];
vert=['constant  ';
      'price     ';
          'sugar     ';
          'mushy     '];

disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	% create weight matrix

invA = inv([IV'*IV]);

disp(['The dimensions of object invA  are ' num2str(size(invA)) ])




disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	% Logit results and save the mean utility as initial values for the
%search below.

	% compute the outside good market share by market:

temp = cumsum(s_jt);
sum1 = temp(cdindex,:);
sum1(2:size(sum1,1),:) = diff(sum1);
outshr = 1.0 - sum1(cdid,:);

disp(['Object temp dimensions: ' num2str(size(temp)) ])
disp(['Object sum1 dimensions: ' num2str(size(sum1)) ])
disp(['Object outshr dimensions: ' num2str(size(outshr)) ])




y = log(s_jt) - log(outshr);
mid = x1'*IV*invA*IV';
t = inv(mid*x1)*mid*y;
mvalold = x1*t;
oldt2 = zeros(size(theta2));
mvalold = exp(mvalold);

 	%s_jt and outshr and x1 are  data that Nevo provides

disp(['Object y dimensions: ' num2str(size(y)) ])
disp(['Object mid dimensions: ' num2str(size(mid)) ])
disp(['Object t dimensions: ' num2str(size(t)) ])
disp(['Object mvalold dimensions: ' num2str(size(mvalold)) ])
disp(['Object oldt2 dimensions: ' num2str(size(oldt2)) ])


	%the next command creates a new file, mvalold.mat, with mvaoldold and oldt2 in it, and then clears out the old mvalold from memory. 

save mvalold mvalold oldt2
clear mid y outshr t oldt2 mvalold temp sum1

vfull = v(cdid,:);
dfull = demogr(cdid,:);

	% the matrix v has 80 iid normal random numbers for each of the 94 observations
 %the matrix demogr has random draws from the CPS for 20 individuals per obs.

disp(['Object vfull dimensions: ' num2str(size(vfull)) ])
disp(['Object dfull dimensions: ' num2str(size(dfull)) ])


disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	%options = foptions; obsolete matlab
	%options(2) = 0.1;obsolete matlab
	%options(3) = 0.01;obsolete matlab
	%foptions is obsolete in Matlab 6. optimset is its replacement.

	%Below,   set the maximum number of iterations for the main optimization command. Note .that the number of iterations for the meanval.m calculation is separate, and unconstrained. 

	%I can't get the Jacobian-Gradobj to work, so I turned it off. It generates a gigantic value for df in gmmobjg.m, and sends the next iteration into outer space. The problem is probably in the gmmobjg.m file, since Bronwyn Hall did not change the jacob.m file. Or, it might be in the meanval.m file somehow. 

options = optimset('GradObj','off','MaxIter',150,'Display','iter','TolFun',0.1,'TolX',0.01);


tic

	% tic starts the clock to measure  the time it takes for optimization

	%Use either (1) or (2) below, but not both, to do the minimization. Use % to comment out the command you are not using. 

	% (1) The following line computes the estimates using a Quasi-Newton method  with an   *analytic* gradient   The old Matlab command fminu is now fminunc, in versions 6 and  7. What is in  Nevo is the obsolete command:

  % OBSOLETE: % [theta2,  options] = fminu('gmmobj',theta2,options,'gradobj')

         %   The new command is below.    Fminunc is  from the  "optimization toolbox", not in the basic  matlab. What you  should use now for the Quasi-Newton method is the next line: 

  %      [theta2,fval,exitflag,output] = fminunc('gmmobjg',theta2,options);

        %gmmobjg.m  is a  separate matlab file created by Bronwyn Hall to %replace  Nevo's gmmobj.m and gradobj.m with one file for Matlab 6 or 7.
 %gmmobjg calls up meanval.m

	 % (2) The following command computes the estimates using a simplex search method.
	%The old Matlab command fmins is now fminsearch, in version 7.

    [theta2,fval,exitflag,output]  = fminsearch('gmmobjg',theta2, options)


	%meanval.m has the command, "disp(['# of iterations for delta convergence:  ' num2str(i)])", which shows the iterations it takes to converge in computing the mean utility for each iteration of the main optimization routine. 

	%gmmobjg.m has the command "disp(['GMM objective:  ' num2str(f1)])", which shows the value of the moment expression at each meanval.m iteration

	%Either (1) or (2) shows various things at each optimization iteration

comp_t = toc/60;

	%toc is  the elapsed optimization time that tic started


disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	% computing the s.e.

vcov = var_cov(theta2);
se = sqrt(diag(vcov));

	%var_cov is a separate matlab file

 

disp(['Object vcov dimensions: ' num2str(size(vcov)) ])
disp(['Object se dimensions: ' num2str(size(se)) ])



theta2w = full(sparse(theti,thetj,theta2));
t = size(se,1) - size(theta2,1);
se2w = full(sparse(theti,thetj,se(t+1:size(se,1))));

disp(['Object theta2w dimensions:     ' num2str(size(theta2w)) ])
disp(['Object t dimensions:     ' num2str(size(t)) ])
disp(['Object se2w dimensions:     ' num2str(size(se2w)) ])


disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	% the MD estimates

omega = inv(vcov(2:25,2:25));
xmd = [x2(1:24,1) x2(1:24,3:4)];
ymd = theta1(2:25);

disp(['Object omega dimensions: ' num2str(size(omega)) ])
disp(['Object xmd dimensions: ' num2str(size(xmd)) ])
disp(['Object ymd dimensions: ' num2str(size(ymd)) ])


beta = inv(xmd'*omega*xmd)*xmd'*omega*ymd;
resmd = ymd - xmd*beta;
semd = sqrt(diag(inv(xmd'*omega*xmd)));
mcoef = [beta(1); theta1(1); beta(2:3)];
semcoef = [semd(1); se(1); semd];

disp(['Object omega dimensions: ' num2str(size(omega)) ])
disp(['Object beta dimensions: ' num2str(size(beta)) ])
disp(['Object resmd dimensions: ' num2str(size(resmd)) ])
disp(['Object semd dimensions: ' num2str(size(semd)) ])
disp(['Object mcoef dimensions: ' num2str(size(mcoef)) ])
disp(['Objectsemcoef dimensions: ' num2str(size(semcoef)) ])


Rsq = 1-((resmd-mean(resmd))'*(resmd-mean(resmd)))/((ymd-mean(ymd))'*(ymd-mean(ymd)));
Rsq_G = 1-(resmd'*omega*resmd)/((ymd-mean(ymd))'*omega*(ymd-mean(ymd)));
Chisq = size(id,1)*resmd'*omega*resmd;


disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

diary results.txt
disp(horz)
disp('  ')
for i=1:size(theta2w,1)
     disp(vert(i,:))
     disp([mcoef(i) theta2w(i,:)])
     disp([semcoef(i) se2w(i,:)])
end

disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

	%disp(['GMM objective:  ' num2str(options(8))]) old matlab version, obsolete

disp(['GMM objective:  ' num2str(fval)])
disp(['MD R-squared:  ' num2str(Rsq)])
disp(['MD weighted R-squared:  ' num2str(Rsq_G)])

	%disp(['# of objective function evaluations:  ' num2str(options(10))]) old matlab, obsolete

disp(['run time (minutes):  ' num2str(comp_t)])

 

diary off

disp('                                                                    ')
disp('**********************************************************************')
disp('                                                                   ')

delete gmmresid.mat 
delete mvalold.mat











