function PS3_solution

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

% load data. see description in readme.txt
load ps3
load iv
IV1= [iv(:,2:21) x2(:,1) x2(:,3:end)];
IV2 = [iv(:,2:21) x1(:,2:25)];

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 %

cdid = kron([1:nmkt]',ones(nbrn,1));  % this vector relates each observation to the market it is in%
cdindex = [nbrn:nbrn:nbrn*nmkt]';     % this vector 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. %


          
horz=['OLS no dum     OLS    IV no dum     IV'];
vert=['constant  ';
      'price     ';
          'sugar     ';
          'mushy     '];

% create weight matrix
invA1 = inv([IV1'*IV1]);
invA2 = inv([IV2'*IV2]);


% Question 1: Generate regression results

% 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,:);

y = log(s_jt) - log(outshr);

beta_ols_nodummies = inv(x2'*x2)*x2'*y;
res = y - x2*beta_ols_nodummies;
%vcov_ols_nodummies = 1/size(x2,1)*res'*res*inv(x2'*x2);
vcov_ols_nodummies = inv(x2'*x2)*(x2'.*(ones(size(x2,2),1)*res')*(res*ones(size(x2,2),1)'.*x2))*inv(x2'*x2);
se_ols_nodummies = sqrt(diag(vcov_ols_nodummies));

beta_ols = inv(x1'*x1)*x1'*y;
res = y - x1*beta_ols;
%vcov_ols = 1/size(x1,1)*res'*res*inv(x1'*x1);
vcov_ols = inv(x1'*x1)*(x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1))*inv(x1'*x1);

% the Min-Distance estimates for sugar and mushy
omega = inv(vcov_ols(2:25,2:25));
xmd = [x2(1:24,1) x2(1:24,3:4)];
ymd = beta_ols(2:25);

beta_temp = inv(xmd'*omega*xmd)*xmd'*omega*ymd;
resmd = ymd - xmd*beta_temp;
%se_temp = sqrt(diag(1/21*resmd'*resmd*inv(xmd'*omega*xmd)));
se_temp = sqrt(diag(inv(xmd'*omega*xmd)*(xmd'.*(ones(size(xmd,2),1)*resmd')*omega*(resmd*ones(size(xmd,2),1)'.*xmd))*inv(xmd'*omega*xmd)));
beta_ols_md = [beta_temp(1); beta_ols(1); beta_temp(2:3)];
se_ols_md = [se_temp(1); sqrt(vcov_ols(1,1)); se_temp(2:3)];

% IV regressions
beta_IV_nodummies = inv(x2'*IV1*invA1*IV1'*x2)*x2'*IV1*invA1*IV1'*y;
res = y - x2*beta_ols_nodummies;
IVres = IV1.*(res*ones(1,size(IV1,2)));
b = IVres'*IVres;
vcov_IV_nodummies = inv(x2'*IV1*invA1*IV1'*x2)*x2'*IV1*invA1*b*invA1*IV1'*x2*inv(x2'*IV1*invA1*IV1'*x2);
se_IV_nodummies = sqrt(diag(vcov_IV_nodummies));

beta_IV = inv(x1'*IV2*invA2*IV2'*x1)*x1'*IV2*invA2*IV2'*y;
res = y - x2*beta_ols_nodummies;
IVres = IV2.*(res*ones(1,size(IV2,2)));
b = IVres'*IVres;
vcov_IV = inv(x1'*IV2*invA2*IV2'*x1)*x1'*IV2*invA2*b*invA2*IV2'*x1*inv(x1'*IV2*invA2*IV2'*x1);

% the Min-Distance estimates for sugar and mushy
omega = inv(vcov_IV(2:25,2:25));
xmd = [x2(1:24,1) x2(1:24,3:4)];
ymd = beta_IV(2:25);

beta_temp = inv(xmd'*omega*xmd)*xmd'*omega*ymd;
resmd = ymd - xmd*beta_temp;
%se_temp = sqrt(diag(1/21*resmd'*resmd*inv(xmd'*omega*xmd)));
se_temp = sqrt(diag(inv(xmd'*omega*xmd)*(xmd'.*(ones(size(xmd,2),1)*resmd')*omega*(resmd*ones(size(xmd,2),1)'.*xmd))*inv(xmd'*omega*xmd)));
beta_IV_md = [beta_temp(1); beta_IV(1); beta_temp(2:3)];
se_IV_md = [se_temp(1); sqrt(vcov_IV(1,1)); se_temp(2:3)];

[beta_ols_nodummies se_ols_nodummies beta_ols_md se_ols_md ...
    beta_IV_nodummies se_IV_nodummies beta_IV_md se_IV_md]

fid = fopen('results.txt','wt');
fprintf(fid,'Question 1\n');
fprintf(fid,'------------\n');
fprintf(fid,'  \n');
fprintf(fid,horz);
fprintf(fid,'   \n');
disp('  ')
for i=1:4
     fprintf(fid, vert(i,:));
     fprintf(fid,'   \n');
     fprintf(fid, '%8.3f %10.3f %8.3f %12.3f\n',[beta_ols_nodummies(i) beta_ols_md(i) beta_IV_nodummies(i) beta_IV_md(i)])
     fprintf(fid, '%8.3f %10.3f %8.3f %12.3f\n',[se_ols_nodummies(i) se_ols_md(i) se_IV_nodummies(i) se_IV_md(i)])

end
fprintf(fid,'  \n')
fprintf(fid,'  \n')

alpha = beta_IV_md(2);
brand = floor(id/100000);
company = floor(brand/1000);
pre_merger_dum = cr_dum(company);
p = x1(:,1);
markup_single = zeros(24*94,1);
markup_multi = zeros(24*94,1);
markup_joint = zeros(24*94,1);
elas = zeros(24*94,24);
omega = zeros(24*94,24);

% compute demand elasticities

for t = 1:nmkt
    shares = s_jt((t-1)*24 + 1:t*24);
    own_dum = pre_merger_dum((t-1)*24 + 1:t*24,:);
    o1 = alpha*shares*shares';
    o2 = diag(alpha*shares);
    omega((t-1)*24 + 1:t*24,:) = (o2-o1);
    elas((t-1)*24 + 1:t*24,:) = omega((t-1)*24 + 1:t*24,:).*(1./shares*p((t-1)*24 + 1:t*24,1)');
    markup_single((t-1)*24 + 1:t*24,1) = -inv(omega((t-1)*24 + 1:t*24,:).*eye(24))*shares;
    markup_multi((t-1)*24 + 1:t*24,1) = -inv(omega((t-1)*24 + 1:t*24,:).*(own_dum*own_dum'))*shares;
    markup_joint((t-1)*24 + 1:t*24,1) = -inv(omega((t-1)*24 + 1:t*24,:).*ones(24))*shares;
end


fprintf(fid,'Question 2\n')
fprintf(fid,'------------\n')
fprintf(fid,'Demand elasticities\n')
fprintf(fid,' \n')
fprintf(fid,' \n')
for i = 1:24 
    fprintf(fid,'%4.3f  %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f\n ',[median(elas(i:24:end,:))]);
end

fclose(fid);

% compute marginal costs
mc_hat_single = p - markup_single;
mc_hat_multi = p - markup_multi;
mc_hat_joint = p - markup_joint;

% compute margins
margin_single = (p - mc_hat_single)./p*100;
margin_multi = (p - mc_hat_multi)./p*100;
margin_joint = (p - mc_hat_joint)./p*100;

diary results.txt
disp([' '])
disp('Question 3')
disp('------------')
disp('markups and margins')
disp([' '])
disp(['              Single product     Multi product       Joint'])
disp(['    brand     markup   margin   markup    margin   markup    margin'])
disp(['---------------------------------------------------------------------'])
disp([brand(1:24)/1000 median(reshape(markup_single,24,94)')' mean(reshape(margin_single,24,94)')' ...
                 median(reshape(markup_multi,24,94)')' mean(reshape(margin_single,24,94)')' ...
                 median(reshape(markup_joint,24,94)')' mean(reshape(margin_joint,24,94)')'])

disp([' '])
disp('Question 4')
disp('------------')
disp('marginal costs')
disp([' '])
disp(['    brand     Single   Multi    Joint'])
disp(['----------------------------------------'])
disp([brand(1:24)/1000 median(reshape(mc_hat_single,24,94)')'  ...
                 median(reshape(mc_hat_multi,24,94)')' ...
                 median(reshape(mc_hat_joint,24,94)')' ])
             
diary off

% merger of Post & Nabisco

% new ownership matrix
company_post = company.*(1-(company==6))+(company==6)*3;
post_dum = cr_dum(company_post);

% solve for post merger price equilbrium
expall = exp(y - alpha*p);
first = 1;
for i=1:94
   last = i*24;
   p_star(first:last,1) = fsolve(@bert_eq, p(first:last), optimset('Display','off'), expall(first:last), alpha, mc_hat_multi(first:last), post_dum(first:last,:));
   first = last + 1;
end   

diary results.txt
disp([' '])
disp('Question 5')
disp('------------')
disp('')
disp ('Post & Nabisco')
disp ('--------------')
disp('')
disp (['    brand     price     p_star   %change '])
disp(['------------------------------------------'])
disp([brand(1:24)/1000 median(reshape(full(p),24,94)')' median(reshape(full(p_star),24,94)')' median(reshape(p_star./p*100-100,24,94)')']) 
diary off


% merger of GM & Quaker

% new ownership matrix
company_post1 = company.*(1-(company==4))+(company==4)*2;
post_dum1 = cr_dum(company_post1);

% solve for post merger price equilbrium
first = 1;
for i=1:94
   last = i*24;
   p_star1(first:last,1) = fsolve(@bert_eq, p(first:last), optimset('Display','off'), expall(first:last), alpha, mc_hat_multi(first:last), post_dum1(first:last,:));
   first = last + 1;
end   

diary results.txt
disp('')
disp ('GM & Quaker')
disp ('--------------')
disp('')
disp (['    brand     price     p_star   %change '])
disp(['------------------------------------------'])
disp([brand(1:24)/1000 median(reshape(full(p),24,94)')' median(reshape(full(p_star1),24,94)')' median(reshape(p_star1./p*100-100,24,94)')']) 
diary off

function f=cr_dum(long_id)

% cr_dum	This function creates a set of dummies for each of the values defined by long_id

b = sort(long_id);
b1 = [1;diff(b)];
b2 = b(b1>0);
clear b1 b
f = sparse(zeros(size(long_id,1), size(b2,1)));
for i = 1:size(b2,1)
	f(:,i) = sparse((long_id==b2(i)));
end

function f = bert_eq(p, expall, alpha, mc_hat, own_dum)

eg = expall.*exp(alpha*p);
shares = eg./(ones(size(eg,1),1)*(1+sum(eg)));
o1 = alpha*shares*shares';
o2 = diag(alpha*shares);
omega = (o2-o1).*(own_dum*own_dum');

f = omega*(p - mc_hat) + shares ;