%%program on multiproduct firms
%%given parameters 
clear
close all
alphabarsigma=2.66  %alpha*(sigma-1)
sigma=8.49
delta=-1.38
theta=8.28
MaxPercCount=40000
Fixed=80
AvSalesData=3170;   %average sales in the data for the calibration of fixed costs
Fixed=AvSalesData/(theta/(sigma-1)*sigma)*(theta/(sigma-1)-1)
power=(delta)*(1-theta/(sigma-1))-[alphabarsigma]/(sigma-1)*theta
power2=(delta)*(-theta/(sigma-1))-[alphabarsigma]/(sigma-1)*theta
constantsum=0
constantsum2=0
constantsum3=0
for i=1:20000
    constantsum=i^power+constantsum;
    constantsum2=i^power2+constantsum2;
    constantsum3=i^[-[alphabarsigma]/(sigma-1)]+constantsum3;
end
avsales=constantsum*Fixed*sigma/(1-(sigma-1)/theta)
avgoods=constantsum2

%Goods distribution data, USA
%----------------------------
GDistr=[1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	3	3	3	3	3	3	3	3	3	3	4	4	4	4	4	4	5	5	5	6	6	7	8	8	10	12	15	21	39;
0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1	0.11	0.12	0.13	0.14	0.15	0.16	0.17	0.18	0.19	0.2	0.21	0.22	0.23	0.24	0.25	0.26	0.27	0.28	0.29	0.3	0.31	0.32	0.33	0.34	0.35	0.36	0.37	0.38	0.39	0.4	0.41	0.42	0.43	0.44	0.45	0.46	0.47	0.48	0.49	0.5	0.51	0.52	0.53	0.54	0.55	0.56	0.57	0.58	0.59	0.6	0.61	0.62	0.63	0.64	0.65	0.66	0.67	0.68	0.69	0.7	0.71	0.72	0.73	0.74	0.75	0.76	0.77	0.78	0.79	0.8	0.81	0.82	0.83	0.84	0.85	0.86	0.87	0.88	0.89	0.9	0.91	0.92	0.93	0.94	0.95	0.96	0.97	0.98	0.99];

%mean scope sales percentiles data (mean scope, mean scale, percentiles)
%----------------------------
MeanScScale=[3.494972	3.517196	3.541363	3.563022	3.579054	3.604302	3.622973	3.650976	3.670074	3.694227	3.718559	3.74016	3.765291	3.788669	3.808446	3.834796	3.856371	3.878859	3.905101	3.932746	3.958249	3.990558	4.018711	4.03412	4.056314	4.090791	4.118755	4.152821	4.175676	4.214253	4.253358	4.245301	4.283262	4.29332	4.327273	4.370259	4.404762	4.421513	4.455544	4.494418	4.537297	4.525563	4.561207	4.605233	4.612044	4.664505	4.723123	4.764382	4.804239	4.843611	4.873541	4.920582	4.917568	4.978621	5.004228	5.073487	5.132646	5.20362	5.275676	5.335968	5.418962	5.502078	5.570819	5.679229	5.711712	5.772222	5.871306	5.979371	5.98075	6.087866	6.19027	6.27933	6.362268	6.466987	6.574813	6.760052	6.866216	7.019718	7.064801	7.035494	7.273906	7.549488	7.763964	7.96	8.287449	8.62851	8.990741	9.311721	9.756757	10.33823	11.02589	11.76619	12.32389	12.44907	13.27568	15.10323	16.16129	18.44086	24.29032	32.3871;
906.9677	910.0925	913.1526	917.0062	922.4559	925.6845	930.4384	933.2746	938.5538	942.7075	946.9882	952.14	956.2225	961.2818	967.4455	972.1368	978.2413	984.3212	989.2716	994.4697	1000.436	1004.918	1010.691	1019.923	1027.267	1032.196	1039.045	1044.638	1053.346	1058.393	1063.138	1080.565	1086.704	1100.301	1108.161	1114.086	1122.004	1135.411	1144.818	1153.416	1161.444	1184.06	1194.259	1203.419	1222.906	1230.924	1237.928	1250.115	1262.543	1276.527	1293.719	1307.115	1334.732	1345.031	1366.703	1377.411	1391.811	1404.039	1417.067	1433.288	1445.665	1459.276	1478.016	1487.673	1518.855	1542.918	1559.684	1575.892	1622.469	1642.831	1666.777	1695.03	1729.272	1760.491	1793.988	1810.017	1851.152	1881.263	1947.811	2041.781	2065.535	2085.8	2130.63	2184.823	2215.567	2253.825	2298.485	2367.019	2420.916	2455.873	2496.302	2555.731	2691.625	2974.288	3160.905	3199.786	3565.507	3904.404	4027.632	4948.476;
100	99.02692	98.02141	97.0159	96.01038	95.00487	94.03179	93.02628	92.02076	91.01525	90.00973	89.00422	88.03114	87.02563	86.02011	85.0146	84.00909	83.00357	82.03049	81.02498	80.01946	79.01395	78.00844	77.00292	76.02985	75.02433	74.01881	73.0133	72.00779	71.00227	70.02919	69.02368	68.01817	67.01265	66.00714	65.00163	64.02854	63.02303	62.01752	61.012	60.00649	59.00098	58.0279	57.02238	56.01687	55.01136	54.00584	53.00032	52.02725	51.02173	50.01622	49.0107	48.00519	47.03211	46.0266	45.02108	44.01557	43.01006	42.00454	41.03146	40.02595	39.02044	38.01492	37.00941	36.00389	35.03082	34.0253	33.01979	32.01427	31.00876	30.00324	29.03016	28.02465	27.01914	26.01362	25.00811	24.0026	23.02952	22.024	21.01849	20.01297	19.00746	18.00195	17.02887	16.02335	15.01784	14.01233	13.00681	12.0013	11.02822	10.02271	9.017191	8.011677	7.006163	6.000649	5.027571	4.022057	3.016542	2.011028	1.005514];

MeanSalesData=MeanScScale(1,:).*MeanScScale(2,:);
%%%first calculate the average number of goods above a threshold
%----------------------------
for i=1:MaxPercCount
    A(i)=unifrnd(0,1);
end
Perc=sort(A);


for i=1:MaxPercCount
    
    Good(i)=floor((1-Perc(i))^(-(sigma-1)/theta/(delta+alphabarsigma)));
    
end




%%%average sales across firms
%  construct the summation term (H(g))
%----------------------------
for i=1:MaxPercCount
    Summation(1)=0;
    for j=1:Good(i)
    Summation(j+1)=Summation(j)+(1/j^(alphabarsigma));
    end
    SummationMat(i)=Summation(Good(i)+1);
end


%now construct sales for different sales percentiles
%----------------------------
for i=1:MaxPercCount  
    Sales(i)=sigma*Fixed*(1-Perc(i))^(-(sigma-1)/theta)*SummationMat(i);
end


for i=1:MaxPercCount
    AvGood(i)=mean(Good(i:end));
    AvSales(i)=mean(Sales(i:end));
end
%generate a super big firm (weighted by distribution)
%----------------------------
SalesBig=sigma*Fixed*(1-.0000001)^(-(sigma-1)/theta)*constantsum3

for i=1:MaxPercCount  
    AvSalesAcross(i)=mean(Sales(i:end))/mean(Good(i:end));
    
end

%creating the percentiles and renormalizing all the vectors in terms of
%percentiles
%----------------------------
for i=1:99
    Perc2(i)=fix(i/100*MaxPercCount);
    AvGood2(i)=AvGood(Perc2(i));
    Good2(i)=Good(Perc2(i));
    Sales2(i)=Sales(Perc2(i));
    AvSalesAcross2(i)=AvSalesAcross(Perc2(i));
    AvSales2(i)=AvSales(Perc2(i));
end



%----------------------------
%% Plotting various distributions
%----------------------------



%Plot Mean Scope Scale for Percentiles ranked by sales size
%----------------------------
scatter(1-Perc2/MaxPercCount,AvGood2,'r')

set(gca,'yscale','log','YTick',[1 10 100 1000],...
  'YTickLabel',{'','10','100','1000'},'XTick',[0.01 .1 1],...
  'XTickLabel',{'0.01','.1','1'})
set(gca,'xscale','log')
xlabel({'Firm Percentile'})
ylabel({'Products and Sales Per Product in the USA'})

hold on
scatter((1-Perc2/MaxPercCount),AvSalesAcross2,'r')
scatter((MeanScScale(3,1:100))/100,MeanScScale(1,1:100),'c','Marker','pentagram')
scatter((MeanScScale(3,1:100))/100,MeanScScale(2,1:100),'c','Marker','pentagram')



%Plot Goods Distributions. Model and Data
%----------------------------
pause
close all
Percentile=Perc2/MaxPercCount;
scatter(Percentile,Good2,'r')
set(gca,'yscale','log','YTick',[1 10 100],...
  'YTickLabel',{'1','10','100'})
xlabel({'Firm Percentile'})
ylabel({'Products in the USA'})

 hold on
scatter((GDistr(2,1:99)),GDistr(1,1:99),'c','Marker','pentagram')




%%plot sales Distribution
%----------------------------
    pause 
    close all
    
    scatter(Percentile,Sales2,'r')
     set(gca,'yscale','log','YTick',[1 10 100 1000 10000],...
      'YTickLabel',{'','10','100','1000','10000'},'XTick',[0.01 .1 1],...
  'XTickLabel',{'0.01','.1','1'})
     xlabel({'Firm Percentile'})
    ylabel({'Sales in the USA'})

%%plot sales Distribution
%----------------------------   
     pause 
    close all
    
    scatter(1-Percentile,AvSales2,'r')
     set(gca,'yscale','log','YTick',[100 1000 10000 100000],...
      'YTickLabel',{'100','1000','10000','100000'},'XTick',[0.01 .1 1],...
  'XTickLabel',{'0.01','.1','1'})
set(gca,'xscale','log')
     xlabel({'Firm Percentile'})
    ylabel({'Average Sales in the USA'})
    hold on 
    scatter((MeanScScale(3,1:100))/100,MeanSalesData,'c')
    
    
%----------------------------
%% Empirical Exercises on the Model Generated Data
%----------------------------

    %regression in the mean sales per firm
    onevec=ones(1,MaxPercCount);  %generating a vector of ones to have a constant
    Goods2=[onevec; log(Good)]'; %generating RHS of the regression, the constant and the goods per firm
    MeanFirmSales=Sales./Good;   %mean sales per product for each firm
    DataVector=[log(Sales./Good) ; log(Good)]';  %Putting data from the model in a vector
    regress(log(MeanFirmSales'),Goods2) %running a regression on the goods and a constant

    
%----------------------------    
%% Saving Variables
%----------------------------
   Data=[Percentile; Good2; Sales2; AvGood2; AvSalesAcross2; AvSales2]
   %variables in the Data (as ordered): Percentiles of the firm, Goods for
   %firms of a given percentile, sales of firms of given percentile,
   %Average Number of Goods for Firms of a Percentile and Above, Scope
   %Weighted Firm Scale for firms from a percentile and above
   
   Data2=Data'
    csvwrite('variables.csv',Data2)
   % save('variables2.mat','Percentile')
    %save('variables.mat','Percentile','AvGood2','Good2','AvSalesAcross2','Sales2')