A simulation of Burdett-Judd competition

Lucas Herrenbrueck, Simon Fraser University
Earliest version: 2013. This version: 2017

Every consumer has one unit of money to spend, and is willing to spend all of it as long as the price is below maxp (see for example: Head & Kumar 2005). Consumers receive a random selection of quotes and then buy from the firm with the lowest price.

For the simulation, I use various procedures (“models” in the paper). There is always a finite number of firms. In the first and second procedure, firms receive a signal about potential profits at one randomly picked price; in the first procedure, they update sequentially, but in the second procedure, they update all at once. The third procedure is like the first but there is cost dispersion. These procedures are repeated with menu costs. In the last procedure, firms update (sequentially) by maximizing profits over a discrete grid of possible prices.

ρ =1; (* non-zero *)
η = 2; (* must satisfy ρη>-1 *)
a[F_]:=(η/(1-(((1+ρ η))^(-1/ρ))))(((1+ρ η F))^(-1-1/ρ)) (* sales as a function of rank, if continuous *)
J[F_]:=((1-(((1+η ρ F))^(-1/ρ)))/(1-(((1+η ρ))^(-1/ρ)))) (* cumulative sales *);

S[x_]:=Min[1,Max[0,1+((1-(((((x-x (((1+η ρ))^(-1/ρ)))/η)))^(-(ρ/(1+ρ)))))/(η ρ)) ]](* benchmark sales cdf *)
F[p_]:=Min[1,Max[0,Re[((-1+(((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/(η ρ)) ]]](* benchmark price cdf *)

cost = 1;  (* must be non-zero *)
maxp = 2;

nprices=10;
pgrid = Table[cost +(maxp-cost)(i/nprices),{i,nprices}]; (* this excludes the cost itself, which firms would never want to charge *)

profitsdiscrete[fgrid_(*share of all firms per price*),pold_(*old price*),nfirms_]:=Table[With[{Fgrid=FoldList[Plus,0,fgrid+Table[Boole[i==pnew]-Boole[i==pold],{i,nprices}]/nfirms]},(1-(cost/(pgrid[[pnew]])))((J[Fgrid[[pnew+1]]]-J[Fgrid[[pnew]]])/(fgrid[[pnew]]+Boole[pold≠pnew]/nfirms))],{pnew(*new price*),nprices}]

Benchmark (modified Burdett-Judd, continuum vs discrete number of firms)

Plot[{1-S[x],F[x]},{x,0,3.5},PlotStyle->{{Thick,Dashed,Darker@Blue},{Thick,Darker@Red}},AxesLabel->{(*Style[" Price,\n",Darker@Red]Style["Sales",Darker@Blue]*)None,"Rank"},BaseStyle->{FontFamily->"Helvetica",Medium,Bold},ImageSize->350,Epilog->{Text[Style[" CDF of prices",Darker@Red],{2.5,.9}],Text[Style["complementary CDF of sales",Darker@Blue],{2.5,.25}]}]

(* Implied density of prices *)
Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}]

Profits for sampling N firms with replacement (which enables k->∞)

Transpose@Table[{x,(1-cost/maxp)x (J[1.]-J[1-1/x])},{x,{2,6,16,40,100,625}}]//TableForm

Update rules for Procedures 1-3

Procedure 1: identical costs, random-order updating by binary comparison

Procedure 2: identical costs, wholesale updating by binary comparison based on the old price distribution

Procedure 3: cost dispersion, random-order updating by binary comparison

update1[{
prices_(*price choice per firm*),
j_(*index of the updating firm*)}]:=
With[{draw=RandomReal[{cost,maxp}],oldrank=(Ordering@Ordering[prices])[[j]]},With[{newrank=(Ordering@Ordering@ReplacePart[prices,j->draw])[[j]]},
If[(1-(cost/draw))nfirms(J[(newrank/nfirms)]-J[((newrank-1)/nfirms)])≤(1-(cost/(prices[[j]])))nfirms(J[(oldrank/nfirms)]-J[((oldrank-1)/nfirms)]),
{prices,Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*nothing changes*),
{ReplacePart[prices,j->draw],Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*ruling out j picking twice in a row*)
]]]

update2[{
prices_(*price choice per firm*),
j_(*empty here*)}]:=
With[{draw=RandomReal[{cost,maxp},nfirms],oldranking=Ordering@Ordering[prices]},With[{newranking=Table[(Ordering@Ordering@ReplacePart[prices,i->draw[[i]]])[[i]],{i,1,nfirms}]},
With[{choices=Table[Boole[(1-(cost/(draw[[i]])))nfirms(J[((newranking[[i]])/nfirms)]-J[((newranking[[i]]-1)/nfirms)])>(1-(cost/(prices[[i]])))nfirms(J[((oldranking[[i]])/nfirms)]-J[((oldranking[[i]]-1)/nfirms)])],{i,nfirms}]},
{prices+choices*(draw-prices),1}
]]]

update3[{
prices_(*price choice per firm*),
j_(*index of the updating firm*)}]:=
With[{draw=RandomReal[{costgrid[[j]],maxp}],oldrank=(Ordering@Ordering[prices])[[j]]},With[{newrank=(Ordering@Ordering@ReplacePart[prices,j->draw])[[j]]},
If[(1-((costgrid[[j]])/draw))nfirms(J[(newrank/nfirms)]-J[((newrank-1)/nfirms)])≤(1-((costgrid[[j]])/(prices[[j]])))nfirms(J[(oldrank/nfirms)]-J[((oldrank-1)/nfirms)]),
{prices,Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*nothing changes*),
{ReplacePart[prices,j->draw],Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*ruling out j picking twice in a row*)
]]]

Animation of Procedure 1 or 2 (N=16)

Interactive animation

SeedRandom[1234]
Clear[nfirms,costgrid]
nfirms=6;
Module[{startcontinuous,results},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
results=NestList[update1,startcontinuous,800];
Animate[With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
GraphicsGrid[
{{Plot[{cdf[p],F[p]},{p,1,2},PlotRange->{{.95,2.05},{-.01,1}},PlotLabel->"Benchmark vs simulated price CDF",ExclusionsStyle->Automatic]},
{Plot[(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),{p,1,2},PlotRange->{{.95,2.05},{0,0.5}},PlotStyle->Red,PlotLabel->"Implied profits of a potential entrant (scaled)",ExclusionsStyle->Automatic]}}]],
{ii,100,200,1},AnimationRate->40,DisplayAllSteps->True,AnimationRunning->False]]

Animation pictures for the paper

See “results”, 5 lines down, for the computation procedure. Put in either “update1” or “update2”.

SeedRandom[1234]
Clear[nfirms,startcontinuous,results,ii]
nfirms=6;
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
results=NestList[update1,startcontinuous,800];

ii=608;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=610;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=616;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=618;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=619;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=620;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=621;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=623;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=624;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

ii=625;
With[{cdf=Sum[UnitStep[#-results[[ii,1,j]]],{j,1,nfirms}]/nfirms&},
Plot[{cdf[p],F[p]},{p,1,2.05},PlotRange->{{.95,2.07},{-.01,1}},ExclusionsStyle->Automatic,AxesLabel->{"price",None},PlotStyle->{{Thick,Darker@Blue},{Thick,Dashed,Darker@Blue}},BaseStyle->Bold,ImageSize->250]]
Plot[Piecewise[{{(1-(cost/p))(nfirms+1)(J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]])/(nfirms+1))]-J[(((Ordering@Ordering@Prepend[results[[ii,1]],p])[[1]]-1)/(nfirms+1))]),p<2}},0],{p,1,2.05},PlotRange->{{.95,2.07},{0,0.4}},Ticks->{Automatic,None},ExclusionsStyle->Automatic,AxesLabel->{"price","profit"},PlotStyle->{Thick,Darker@Red},BaseStyle->Bold,ImageSize->250]

Clear[nfirms,startcontinuous,results,ii]

Simulation of Procedure 1

The “best behaved” one. Update opportunities arrive in random sequential order: a firm is shown what its profits would be at another, randomly selected, price, and switches iff they exceed profits at the current price.

N=2 (special graphs using exact benchmark)

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=2;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{0,"(Subscript[σ, ∞])=0"},{(1/16) (4-3 Log[3]),"(Subscript[σ, 2])"},{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])/2"},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),"(Subscript[π, 2])/2"}},None},GridLines->{{{(1/16) (4-3 Log[3]),{Gray,Dotted}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,DotDashed}}},None},PlotRange->{{0,2(1-cost/maxp)(J[1.]-J[1-1/nfirms])},{0,10nfirms}}]
Show[Plot[(0.25/(((1-p))^2)),{p,1.2,2},PlotRange->{{1.,2},{0,3}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=6

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=6;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

Table[SmoothHistogram[Flatten[simulatedprices],{Automatic,i},PlotStyle->{Thick,Darker@Hue[.3,.5]},PlotLabel->i],{i,{"Biweight","Cosine","Epanechnikov","Gaussian","Rectangular","SemiCircle","Triangular","Triweight"}}]

N=16

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=40

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=40;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=100

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=100;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=625

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=625;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

Simulation of Procedure 2

Like Procedure 1, but everyone gets an update opportunity (simultaneously) every period.

N=6

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=6;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=16

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=40

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=40;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=100

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=100;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

ACFs

acflength=Min[nfirms*10,simulationlength/2];
acflength=20;
acfmeanprofits=Table[{k,Correlation[Drop[meanprofits,k],Take[meanprofits,simulationlength-k]]},{k,0,acflength}];
acfvarprofits=Table[{k,Correlation[Drop[dispprofits,k]^2,Take[dispprofits,simulationlength-k]^2]},{k,0,acflength}];

ListPlot[{acfmeanprofits,acfvarprofits},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},PlotRange->All,BaseStyle->{Bold,Medium}(*,PlotLabel->"Autocorrelation functions for the mean (blue)\nand the cross\[Hyphen]section variance (red) of profits"*)(*,Ticks->{Range[0,13*1,1],Automatic}*)]

Result: Using procedure 1, the ACF of the mean hits 0 around 40 (with 40 firms), and around 100 with 100 firms. That’s reassuring. With procedure 4, it takes about 60 (with 40 firms).

Result: Using procedure 2, the ACF is negative for both profits and dispersion at cyclical point. Seems to be a cycle of length 6. Its amplitude increases strongly with N. (Looks like the cycle length increases with menu costs, too -- but the amplitude is unchanged.)

P1, N=100, no menu cost:

P2, N=100, no menu cost:

Animation of Procedure 3 (N=16)

SeedRandom[1234]
Clear[nfirms,costgrid]
nfirms=16;
costgrid=cost(1+.4((Range[nfirms]-1)/(nfirms-1)));
Module[{startcontinuous,results,meanprices},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
results=NestList[update3,startcontinuous,20000];
meanprices=Mean@Transpose[results][[1]];
Animate[
ListPlot[{Transpose[{costgrid,results[[ii,1]]}],Transpose[{costgrid,meanprices}]},Joined->{False,True},PlotStyle->{{Thick,Darker@Blue},{Darker@Red}},AspectRatio->1,PlotRange->{{.95,2.05},{.95,2.05}},PlotLabel->"Prices vs costs",BaseStyle->Bold],
{ii,1000,2000,5},AnimationRate->40,DisplayAllSteps->True,AnimationRunning->False]]

Simulation of Procedure 3 (N=16)

With cost dispersion (fixed uniform grid).

(OverBar[c])=1.1

SeedRandom[1234]
Clear[nfirms,costgrid,simulationburnin,simulationlength,simulatedprices,simulatedranks,correlations,pricemeans,pricesds,priceupper,pricelower,pricemedian,numchanges]
nfirms=16;
costgrid=cost(1+.1((Range[nfirms]-1)/(nfirms-1)));
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update3,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
correlations=N@MapThread[Correlation,{simulatedranks,Table[ii,{simulationlength},{ii,nfirms}]}];
pricemeans=Mean@simulatedprices;
pricesds=Map[RootMeanSquare,Transpose[simulatedprices-Table[pricemeans,{simulationlength}]]];
priceupper=Thread[Quantile[simulatedprices,.95]];
pricelower=Thread[Quantile[simulatedprices,.05]];
pricemedian=Thread[Quantile[simulatedprices,.5]];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["average correlation = ",Mean[correlations]]
Print["SD of the correlation = ",N[StandardDeviation[correlations]]]
(* this being the RANK correlation of prices and costs *)

Graph of range, with means and 2SD intervals:

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,pricemeans+2pricesds]}],Transpose[{costgrid,pricemeans}],Transpose[{costgrid,MapThread[Max,{costgrid,pricemeans-2pricesds}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->{Medium,Bold},ImageSize->250]

Graph of range, with medians and 95-5 percentile intervals:

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,priceupper]}],Transpose[{costgrid,pricemedian}],Transpose[{costgrid,MapThread[Max,{costgrid,pricelower}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->{Medium,Bold},ImageSize->250]

(OverBar[c])=1.4

SeedRandom[1234]
Clear[nfirms,costgrid,simulationburnin,simulationlength,simulatedprices,simulatedranks,correlations,pricemeans,pricesds,priceupper,pricelower,pricemedian,numchanges]
nfirms=16;
costgrid=cost(1+.4((Range[nfirms]-1)/(nfirms-1)));
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update3,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
correlations=N@MapThread[Correlation,{simulatedranks,Table[ii,{simulationlength},{ii,nfirms}]}];
pricemeans=Mean@simulatedprices;
pricesds=Map[RootMeanSquare,Transpose[simulatedprices-Table[pricemeans,{simulationlength}]]];
priceupper=Thread[Quantile[simulatedprices,.95]];
pricelower=Thread[Quantile[simulatedprices,.05]];
pricemedian=Thread[Quantile[simulatedprices,.5]];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["average correlation = ",Mean[correlations]]
Print["SD of the correlation = ",N[StandardDeviation[correlations]]]
(* this being the RANK correlation of prices and costs *)

Graph of range, with means and 2SD intervals:

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,pricemeans+2pricesds]}],Transpose[{costgrid,pricemeans}],Transpose[{costgrid,MapThread[Max,{costgrid,pricemeans-2pricesds}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->{Medium,Bold},ImageSize->250]

Graph of range, with medians and 95-5 percentile intervals:

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,priceupper]}],Transpose[{costgrid,pricemedian}],Transpose[{costgrid,MapThread[Max,{costgrid,pricelower}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->{Medium,Bold},ImageSize->250]

(OverBar[c])=1.9

SeedRandom[1234]
Clear[nfirms,costgrid,simulationburnin,simulationlength,simulatedprices,simulatedranks,correlations,pricemeans,pricesds,priceupper,pricelower,pricemedian,numchanges]
nfirms=16;
costgrid=cost(1+.9((Range[nfirms]-1)/(nfirms-1)));
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update3,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
correlations=N@MapThread[Correlation,{simulatedranks,Table[ii,{simulationlength},{ii,nfirms}]}];
pricemeans=Mean@simulatedprices;
pricesds=Map[RootMeanSquare,Transpose[simulatedprices-Table[pricemeans,{simulationlength}]]];
priceupper=Thread[Quantile[simulatedprices,.95]];
pricelower=Thread[Quantile[simulatedprices,.05]];
pricemedian=Thread[Quantile[simulatedprices,.5]];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["average correlation = ",Mean[correlations]]
Print["SD of the correlation = ",N[StandardDeviation[correlations]]]
(* this being the RANK correlation of prices and costs *)

Graph of range, with means and 2SD intervals:

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,pricemeans+2pricesds]}],Transpose[{costgrid,pricemeans}],Transpose[{costgrid,MapThread[Max,{costgrid,pricemeans-2pricesds}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->{Medium,Bold},ImageSize->250]

Graph of range, with medians and 95-5 percentile intervals:

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,priceupper]}],Transpose[{costgrid,pricemedian}],Transpose[{costgrid,MapThread[Max,{costgrid,pricelower}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->{Medium,Bold},ImageSize->250]

Simulation with menu cost (Model 4)

Do N=16 (because that’s the number I used for P3). Then compare P1, P2, and P3 with c=1.4, for menu costs of 1/60 (10% of π∞), 1/30 (20%), and 1/15 (40%), with a slight modification for P3 (see below).

Update rules

Procedure 1: identical costs, random-order updating by binary comparison

Procedure 2: identical costs, wholesale updating by binary comparison based on the old price distribution

Procedure 3: cost dispersion, random-order updating by binary comparison

update1mc[{
prices_(*price choice per firm*),
j_(*index of the updating firm*)}]:=
With[{draw=RandomReal[{cost,maxp}],oldrank=(Ordering@Ordering[prices])[[j]]},With[{newrank=(Ordering@Ordering@ReplacePart[prices,j->draw])[[j]]},
If[(1-(cost/draw))nfirms(J[(newrank/nfirms)]-J[((newrank-1)/nfirms)])-d≤(1-(cost/(prices[[j]])))nfirms(J[(oldrank/nfirms)]-J[((oldrank-1)/nfirms)]),
{prices,Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*nothing changes*),
{ReplacePart[prices,j->draw],Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*ruling out j picking twice in a row*)
]]]

update2mc[{
prices_(*price choice per firm*),
j_(*empty here*)}]:=
With[{draw=RandomReal[{cost,maxp},nfirms],oldranking=Ordering@Ordering[prices]},With[{newranking=Table[(Ordering@Ordering@ReplacePart[prices,i->draw[[i]]])[[i]],{i,1,nfirms}]},
With[{choices=Table[Boole[(1-(cost/(draw[[i]])))nfirms(J[((newranking[[i]])/nfirms)]-J[((newranking[[i]]-1)/nfirms)])-d>(1-(cost/(prices[[i]])))nfirms(J[((oldranking[[i]])/nfirms)]-J[((oldranking[[i]]-1)/nfirms)])],{i,nfirms}]},
{prices+choices*(draw-prices),1}
]]]

update3mc[{
prices_(*price choice per firm*),
j_(*index of the updating firm*)}]:=
With[{draw=RandomReal[{costgrid[[j]],maxp}],oldrank=(Ordering@Ordering[prices])[[j]]},With[{newrank=(Ordering@Ordering@ReplacePart[prices,j->draw])[[j]]},
If[(1-((costgrid[[j]])/draw))nfirms(J[(newrank/nfirms)]-J[((newrank-1)/nfirms)])-d≤(1-((costgrid[[j]])/(prices[[j]])))nfirms(J[(oldrank/nfirms)]-J[((oldrank-1)/nfirms)]),
{prices,Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*nothing changes*),
{ReplacePart[prices,j->draw],Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*ruling out j picking twice in a row*)
]]]

Simulation runs with d=0.1*π

P1

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
d=0.1N[(1-cost/maxp)a[1]];
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"benchmark"}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

P2

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
d=0.1N[(1-cost/maxp)a[1]];
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"benchmark"}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

P3 ((OverBar[c])=1.4)

Now, assess the menu cost as a proportion of the minimal profit predicted by the continuum-benchmark, i.e., the profit of the highest-cost firm. The cheapest firm will have a higher expected profit than this.

N[(1-1.4cost/maxp)a[1]]

SeedRandom[1234]
Clear[nfirms,costgrid,simulationburnin,simulationlength,simulatedprices,simulatedranks,correlations,pricemeans,pricesds,numchanges]
d=0.1N[(1-1.4cost/maxp)a[1]];
nfirms=16;
costgrid=cost(1+.4((Range[nfirms]-1)/(nfirms-1)));
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update3mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
correlations=N@MapThread[Correlation,{simulatedranks,Table[ii,{simulationlength},{ii,nfirms}]}];
pricemeans=Mean@simulatedprices;
pricesds=Map[RootMeanSquare,Transpose[simulatedprices-Table[pricemeans,{simulationlength}]]];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["average correlation = ",Mean[correlations]]
Print["SD of the correlation = ",N[StandardDeviation[correlations]]]
(* this being the rank correlation of prices and costs *)

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,pricemeans+2pricesds]}],Transpose[{costgrid,pricemeans}],Transpose[{costgrid,MapThread[Max,{costgrid,pricemeans-2pricesds}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->Bold,ImageSize->250]

Simulation runs with d=0.2*π

P1

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
d=0.2N[(1-cost/maxp)a[1]];
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"benchmark"}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

P2

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
d=0.2N[(1-cost/maxp)a[1]];
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"benchmark"}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

P3 ((OverBar[c])=1.4)

Now, assess the menu cost as a proportion of the minimal profit predicted by the continuum-benchmark, i.e., the profit of the highest-cost firm. The cheapest firm will have a higher expected profit than this.

N[(1-1.4cost/maxp)a[1]]

SeedRandom[1234]
Clear[nfirms,costgrid,simulationburnin,simulationlength,simulatedprices,simulatedranks,correlations,pricemeans,pricesds,numchanges]
d=0.2N[(1-1.4cost/maxp)a[1]];
nfirms=16;
costgrid=cost(1+.4((Range[nfirms]-1)/(nfirms-1)));
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update3mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
correlations=N@MapThread[Correlation,{simulatedranks,Table[ii,{simulationlength},{ii,nfirms}]}];
pricemeans=Mean@simulatedprices;
pricesds=Map[RootMeanSquare,Transpose[simulatedprices-Table[pricemeans,{simulationlength}]]];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["average correlation = ",Mean[correlations]]
Print["SD of the correlation = ",N[StandardDeviation[correlations]]]
(* this being the RANK correlation of prices and costs *)

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,pricemeans+2pricesds]}],Transpose[{costgrid,pricemeans}],Transpose[{costgrid,MapThread[Max,{costgrid,pricemeans-2pricesds}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->Bold,ImageSize->250]

Simulation runs with d=0.4*π

P1

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
d=0.4N[(1-cost/maxp)a[1]];
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update1mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"benchmark"}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

P2

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
d=0.4N[(1-cost/maxp)a[1]];
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update2mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];
numchanges=Count[Flatten@MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices},2],False];

Print["update frequency = ",numchanges/(simulationlength-1.)/nfirms]
Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2.5(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"benchmark"}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,30nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

P3 ((OverBar[c])=1.4)

Now, assess the menu cost as a proportion of the minimal profit predicted by the continuum-benchmark, i.e., the profit of the highest-cost firm. The cheapest firm will have a higher expected profit than this.

N[(1-1.4cost/maxp)a[1]]

SeedRandom[1234]
Clear[nfirms,costgrid,simulationburnin,simulationlength,simulatedprices,simulatedranks,correlations,pricemeans,pricesds,numchanges]
d=0.4N[(1-1.4cost/maxp)a[1]];
nfirms=16;
costgrid=cost(1+.4((Range[nfirms]-1)/(nfirms-1)));
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[update3mc,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
correlations=N@MapThread[Correlation,{simulatedranks,Table[ii,{simulationlength},{ii,nfirms}]}];
pricemeans=Mean@simulatedprices;
pricesds=Map[RootMeanSquare,Transpose[simulatedprices-Table[pricemeans,{simulationlength}]]];
numchanges=Count[MapThread[Equal,{Most@simulatedprices,Rest@simulatedprices}],False];

Print["update frequency = ",numchanges/(simulationlength-1.)]
Print["average correlation = ",Mean[correlations]]
Print["SD of the correlation = ",N[StandardDeviation[correlations]]]
(* this being the rank correlation of prices and costs *)

ListPlot[{Transpose[{costgrid,Map[Min[maxp,#]&,pricemeans+2pricesds]}],Transpose[{costgrid,pricemeans}],Transpose[{costgrid,MapThread[Max,{costgrid,pricemeans-2pricesds}]}]},Joined->False,PlotStyle->{{Darker@Red},{Thick,Darker@Blue}},AxesLabel->{"cost","price"},Filling->{1->{2},3->{2}},FillingStyle->Directive[Opacity[.2],Darker@Blue],AspectRatio->1,PlotRange->{{1,2},{1,2}},PlotRangePadding->.05,BaseStyle->Bold,ImageSize->250]

Simulation with profitmax on a grid (Models A.1 and A.2)

Update rules for Procedures A.1-2

updateA1[{
prices_(*price choice per firm, given as INDEX*),
fgrid_(*share of all firms per price*),
j_(*index of the updating firm*)}]:=
With[{choice=Ordering[profitsdiscrete[fgrid,prices[[j]],nfirms],-1][[1]]},
If[choice==prices[[j]],{prices,fgrid,Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*nothing changes*),
{ReplacePart[prices,j->choice],
ReplacePart[fgrid,{prices[[j]]->fgrid[[prices[[j]]]]-1/nfirms,choice->fgrid[[choice]]+1/nfirms}],
Mod[j+RandomInteger[nfirms-2],nfirms]+1 (*ruling out j again as there wouldn't be anything to update*)
}]]

updateA2[{
prices_(*price choice per firm, given as INDEX*),
fgrid_(*share of all firms per price*),
j_(*index of the updating firm*)}]:=
With[{choice=Ordering[profitsdiscrete[fgrid,prices[[j]],nfirms],-1][[1]]},
If[choice==prices[[j]],{prices,fgrid,Mod[j,nfirms]+1}(*nothing changes*),
{ReplacePart[prices,j->choice],
ReplacePart[fgrid,{prices[[j]]->fgrid[[prices[[j]]]]-1/nfirms,choice->fgrid[[choice]]+1/nfirms}],
Mod[j,nfirms]+1
}]]

Simulation of Procedure A.1

Profitmax on a discrete grid. Updating in random sequential order.

N=16

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits]
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simraw=Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
Transpose[NestList[updateA1,startdiscrete,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
];
simulatedprices=Map[pgrid[[#]]&,simraw[[1]],{2}];
simulatedpmf=simraw[[2]];
simulatedcdf=Map[FoldList[Plus,0,#]&,simulatedpmf];
meanprofits=1./nfirms Table[Sum[(1-(cost/(pgrid[[j]])))(J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]]),{j,nprices}],{i,simulationlength}];
dispprofits=Table[Sqrt[Sum[If[simulatedpmf[[i,j]]==0,0.,simulatedpmf[[i,j]]((((1-(cost/(pgrid[[j]])))(((J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]])/(nfirms simulatedpmf[[i,j]])))-meanprofits[[i]]))^2)],{j,nprices}]],{i,simulationlength}];
(*simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];*)

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=40

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits]
nfirms=40;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simraw=Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
Transpose[NestList[updateA1,startdiscrete,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
];
simulatedprices=Map[pgrid[[#]]&,simraw[[1]],{2}];
simulatedpmf=simraw[[2]];
simulatedcdf=Map[FoldList[Plus,0,#]&,simulatedpmf];
meanprofits=1./nfirms Table[Sum[(1-(cost/(pgrid[[j]])))(J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]]),{j,nprices}],{i,simulationlength}];
dispprofits=Table[Sqrt[Sum[If[simulatedpmf[[i,j]]==0,0.,simulatedpmf[[i,j]]((((1-(cost/(pgrid[[j]])))(((J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]])/(nfirms simulatedpmf[[i,j]])))-meanprofits[[i]]))^2)],{j,nprices}]],{i,simulationlength}];
(*simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];*)

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=100

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits]
nfirms=100;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simraw=Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
Transpose[NestList[updateA1,startdiscrete,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
];
simulatedprices=Map[pgrid[[#]]&,simraw[[1]],{2}];
simulatedpmf=simraw[[2]];
simulatedcdf=Map[FoldList[Plus,0,#]&,simulatedpmf];
meanprofits=1./nfirms Table[Sum[(1-(cost/(pgrid[[j]])))(J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]]),{j,nprices}],{i,simulationlength}];
dispprofits=Table[Sqrt[Sum[If[simulatedpmf[[i,j]]==0,0.,simulatedpmf[[i,j]]((((1-(cost/(pgrid[[j]])))(((J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]])/(nfirms simulatedpmf[[i,j]])))-meanprofits[[i]]))^2)],{j,nprices}]],{i,simulationlength}];
(*simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];*)

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

Simulation of Procedure A.2

Profitmax on a discrete grid. Updating in fixed sequential order.

N=16

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits]
nfirms=16;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simraw=Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
Transpose[NestList[updateA2,startdiscrete,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
];
simulatedprices=Map[pgrid[[#]]&,simraw[[1]],{2}];
simulatedpmf=simraw[[2]];
simulatedcdf=Map[FoldList[Plus,0,#]&,simulatedpmf];
meanprofits=1./nfirms Table[Sum[(1-(cost/(pgrid[[j]])))(J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]]),{j,nprices}],{i,simulationlength}];
dispprofits=Table[Sqrt[Sum[If[simulatedpmf[[i,j]]==0,0.,simulatedpmf[[i,j]]((((1-(cost/(pgrid[[j]])))(((J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]])/(nfirms simulatedpmf[[i,j]])))-meanprofits[[i]]))^2)],{j,nprices}]],{i,simulationlength}];
(*simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];*)

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

(* CYCLE LENGTH *)
Length[meanprofits]/With[{list=Take[Abs@Fourier[meanprofits],{2,1000}]},Ordering[list,-1]]//N
%/(2nfirms)

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=40

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits]
nfirms=40;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simraw=Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
Transpose[NestList[updateA2,startdiscrete,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
];
simulatedprices=Map[pgrid[[#]]&,simraw[[1]],{2}];
simulatedpmf=simraw[[2]];
simulatedcdf=Map[FoldList[Plus,0,#]&,simulatedpmf];
meanprofits=1./nfirms Table[Sum[(1-(cost/(pgrid[[j]])))(J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]]),{j,nprices}],{i,simulationlength}];
dispprofits=Table[Sqrt[Sum[If[simulatedpmf[[i,j]]==0,0.,simulatedpmf[[i,j]]((((1-(cost/(pgrid[[j]])))(((J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]])/(nfirms simulatedpmf[[i,j]])))-meanprofits[[i]]))^2)],{j,nprices}]],{i,simulationlength}];
(*simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];*)

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

(* CYCLE LENGTH *)
Length[meanprofits]/With[{list=Take[Abs@Fourier[meanprofits],{2,1000}]},Ordering[list,-1]]//N
%/(2nfirms)

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

N=100

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits]
nfirms=100;
simulationburnin=1000;
simulationlength=20000; (* must be even *)
simraw=Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
Transpose[NestList[updateA2,startdiscrete,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
];
simulatedprices=Map[pgrid[[#]]&,simraw[[1]],{2}];
simulatedpmf=simraw[[2]];
simulatedcdf=Map[FoldList[Plus,0,#]&,simulatedpmf];
meanprofits=1./nfirms Table[Sum[(1-(cost/(pgrid[[j]])))(J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]]),{j,nprices}],{i,simulationlength}];
dispprofits=Table[Sqrt[Sum[If[simulatedpmf[[i,j]]==0,0.,simulatedpmf[[i,j]]((((1-(cost/(pgrid[[j]])))(((J[simulatedcdf[[i,j+1]]]-J[simulatedcdf[[i,j]]])/(nfirms simulatedpmf[[i,j]])))-meanprofits[[i]]))^2)],{j,nprices}]],{i,simulationlength}];
(*simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];*)

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

(* CYCLE LENGTH *)
Length[meanprofits]/With[{list=Take[Abs@Fourier[meanprofits],{2,1000}]},Ordering[list,-1]]//N
%/(2nfirms)

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]
Show[Plot[-((cost (((((maxp (cost-p) (((1+η ρ))^(1+(1/ρ))))/((cost-maxp) p))))^(ρ/(1+ρ))))/((cost-p) p η (1+ρ))),{p,((cost maxp (((1+η ρ))^(1+(1/ρ))))/(cost+maxp (-1+(((1+η ρ))^(1+(1/ρ)))))),maxp},PlotRange->{{cost,maxp},{0,Automatic}},Filling->Axis,FillingStyle->Directive[Darker@Hue[.3,.5],Opacity[0.2]],PlotStyle->{Dashed,Darker@Hue[.3,.5]}],
SmoothHistogram[Flatten[simulatedprices],PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{cost,maxp},Automatic},AxesOrigin->{1,0}]

Animation of Procedure A.2 (N=40)

nfirms=40;
Module[{startdiscrete},
startdiscrete = With[{prices=RandomInteger[{1,nprices},nfirms]},{prices,Drop[BinCounts[prices],1]/nfirms,1}];
With[{results=NestList[updateA2,startdiscrete,400]},
Animate[
GraphicsGrid[{{
Show[
Plot[F[p],{p,1,2},PlotRange->{{.95,2.05},{-.01,1.01}},PlotLabel->"Benchmark vs simulated price CDF"],
ListPlot[Transpose[{Prepend[pgrid,cost],Prepend[Accumulate[results[[ii,2]]],0]}],PlotStyle->PointSize[.015]]]},{
Show[
ListPlot[Transpose[{Prepend[pgrid,cost],Prepend[results[[ii,2]],0]}],Filling->Axis,Joined->False,PlotRange->{{.95,2.05},{-.01,0.5}},PlotStyle->PointSize[.015],PlotLabel->"Simulated price PMF (blue) and\nimplied profits of a potential entrant (red, scaled)"],
ListPlot[Transpose[{Prepend[pgrid,cost],(*Prepend[Table[With[{Fgrid=Prepend[Accumulate[results[[ii,2]]],0]},(1-(cost/(pgrid[[pnew]])))((J[Fgrid[[pnew+1]]]-J[Fgrid[[pnew]]])/(results[[ii,2,pnew]]+1/nfirms))],{pnew(*new price*),nprices}],0]*)
Prepend[profitsdiscrete[results[[ii,2]],0,nfirms],0.01]}],Joined->False,PlotStyle->{Red,PointSize[.015]}]]
}}],{ii,101,300,1},AnimationRate->80,DisplayAllSteps->True,AnimationRunning->False]]
]

For comparison: randomizing rather than updating

Every firm picks from the continuum-benchmark price distribution every period. Compute dispersion and volatility.

Pricing:

Finv[z_]:=(cost/(1+(((cost-maxp) (((1+η ρ))^(-((1+ρ)/ρ))) (((1+z η ρ))^(1+(1/ρ))))/maxp)))
updateR[{
prices_(*price choice per firm*),
j_(*empty here*)}]:=
With[{draw=Finv/@RandomReal[{0,1},nfirms]},{draw,1}]

N=6:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=6;
simulationburnin=1;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[updateR,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]

N=16:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=16;
simulationburnin=1;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[updateR,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]

N=40:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=40;
simulationburnin=1;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[updateR,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]

N=100:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=100;
simulationburnin=1;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[updateR,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]

N=625:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedprices,simulatedranks,simulatedprofits,meanprofits,dispprofits,numchanges]
nfirms=625;
simulationburnin=1;
simulationlength=20000; (* must be even *)
simulatedprices=Module[{startcontinuous},
startcontinuous=With[{prices=RandomReal[{cost,maxp},nfirms]},{prices,1}];
Transpose[NestList[updateR,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedranks=Map[Ordering@Ordering[#]&,simulatedprices];
simulatedprofits=MapThread[(1-(cost/(#1)))(J[((#2)/nfirms)]-J[((#2-1)/nfirms)])&,{simulatedprices,simulatedranks},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[(1-cost/maxp)a[1]]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(1-cost/maxp)a[1]/nfirms}]
SmoothHistogram[{meanprofits,dispprofits},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(1-cost/maxp)a[1]/nfirms,"(Subscript[π, ∞])"}},None},GridLines->{{{(1-cost/maxp)a[1]/nfirms,{Gray,DotDashed}},{(1-cost/maxp)(J[1.]-J[1-1/nfirms]),{Gray,Dotted}}},None},PlotRange->{{0,2(1-cost/maxp)a[1]/nfirms},{0,100nfirms}}]

For comparison: simulation of Cournot competition

Demand is P = maxp - slope*Q. For comparison with the BJ model above, keep cost and maxp the same (set cost=1 and maxp=2). In contrast to BJ however, firms now compete in quantities not prices. But this is not an essential distinction, and I think (say) simulating quantal-response price competition would look similar to Cournot -- the key distinction is in between an equilibrium of endogenous price dispersion, and a setting where every firm has a unique best response that varies smoothly with other firms’ actions.

Update rules

Procedure 1: identical costs, random-order updating by binary comparison

Procedure 2: identical costs, wholesale updating by binary comparison based on the old price distribution

slope=1;
maxq=maxp/slope;
update1cour[{
quantities_(*quantity choice per firm*),
j_(*index of the updating firm*)}]:=
Module[{draw,oldprice,newprice},
(*draw=RandomReal[{0,maxq-Total[quantities]+quantities[[j]]}];(*maximal quantity that would drive everyone's profits to zero, assuming competitors' quantities are fixed*)*)
draw=RandomReal[{0,maxq/(nfirms+1)*2}];(*alternative: twice the CNE quantity, which is also the quantity that would drive profits to zero IN THE CNE*)
oldprice=maxp-slope Total[quantities];
newprice=maxp-slope Total@ReplacePart[quantities,j->draw];
If[(newprice-cost)draw≤(oldprice-cost)quantities[[j]],
{quantities,Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*nothing changes*),
{ReplacePart[quantities,j->draw],Mod[j+RandomInteger[nfirms-2],nfirms]+1}(*ruling out j picking twice in a row*)
]]
update2cour[{
quantities_(*quantity choice per firm*),
j_(*empty here*)}]:=
With[{draw=RandomReal[{0,maxq/(nfirms+1)*2},nfirms],oldprice=maxp-slope Total[quantities]},With[{newprice=maxp-slope Table[Total@ReplacePart[quantities,i->draw[[i]]],{i,1,nfirms}]},
With[{choices=Table[Boole[(newprice[[i]]-cost)draw[[i]]>(oldprice-cost)quantities[[i]]],{i,nfirms}]},
{quantities+choices*(draw-quantities),1}
]]]

Simulation runs

Simulate Procedure 1 for N=6:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedquantities,simulatedtotals,simulatedprofits,meanprofits,dispprofits]
nfirms=6;
simulationburnin=nfirms^2;
simulationlength=20000; (* must be even *)
simulatedquantities=Module[{startcontinuous},
startcontinuous=With[{quantities=RandomReal[{0,maxq/nfirms},nfirms]},{quantities,1}];
Transpose[NestList[update1cour,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedtotals=Transpose@Table[Total[simulatedquantities,{2}],{nfirms}];
simulatedprofits=MapThread[(maxp-slope#2-cost)#1&,{simulatedquantities,simulatedtotals},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[nfirms(((maxp-cost)/(nfirms+1)))^2/slope]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(((maxp-cost)/(nfirms+1)))^2/slope}]
SmoothHistogram[{meanprofits,dispprofits},{"Standardized", 1/4},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(((maxp-cost)/(nfirms+1)))^2/slope,"(Subscript[π, N])"}},None},PlotRange->{{0,2(((maxp-cost)/(nfirms+1)))^2/slope},{0,100nfirms}}]
Show[SmoothHistogram[Flatten[simulatedquantities],{"Standardized", 1/4},PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{0,maxq/nfirms},Automatic},AxesOrigin->{0,0}]

Table[SmoothHistogram[Flatten[simulatedquantities],{{"Standardized", 1/2},i},PlotStyle->{Thick,Darker@Hue[.3,.5]},PlotLabel->i],{i,{"Biweight","Cosine","Epanechnikov","Gaussian","Rectangular","SemiCircle","Triangular","Triweight"}}]

Simulate Procedure 1 for N=16:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedquantities,simulatedtotals,simulatedprofits,meanprofits,dispprofits]
nfirms=16;
simulationburnin=nfirms^2;
simulationlength=20000; (* must be even *)
simulatedquantities=Module[{startcontinuous},
startcontinuous=With[{quantities=RandomReal[{0,maxq/nfirms},nfirms]},{quantities,1}];
Transpose[NestList[update1cour,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedtotals=Transpose@Table[Total[simulatedquantities,{2}],{nfirms}];
simulatedprofits=MapThread[(maxp-slope#2-cost)#1&,{simulatedquantities,simulatedtotals},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[nfirms(((maxp-cost)/(nfirms+1)))^2/slope]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,4000],Take[dispprofits,4000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(((maxp-cost)/(nfirms+1)))^2/slope}]
SmoothHistogram[{meanprofits,dispprofits},{"Standardized", 1/4},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(((maxp-cost)/(nfirms+1)))^2/slope,"(Subscript[π, N])"}},None},PlotRange->{{0,2(((maxp-cost)/(nfirms+1)))^2/slope},{0,100nfirms}}]
Show[SmoothHistogram[Flatten[simulatedquantities],{"Standardized", 1/4},PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{0,maxq/nfirms},Automatic},AxesOrigin->{0,0}]

Simulate Procedure 2 for N=6:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedquantities,simulatedtotals,simulatedprofits,meanprofits,dispprofits]
nfirms=6;
simulationburnin=nfirms^2;
simulationlength=20000; (* must be even *)
simulatedquantities=Module[{startcontinuous},
startcontinuous=With[{quantities=RandomReal[{0,maxq/nfirms},nfirms]},{quantities,1}];
Transpose[NestList[update2cour,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedtotals=Transpose@Table[Total[simulatedquantities,{2}],{nfirms}];
simulatedprofits=MapThread[(maxp-slope#2-cost)#1&,{simulatedquantities,simulatedtotals},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[nfirms(((maxp-cost)/(nfirms+1)))^2/slope]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(((maxp-cost)/(nfirms+1)))^2/slope}]
SmoothHistogram[{meanprofits,dispprofits},{"Standardized", 1/4},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(((maxp-cost)/(nfirms+1)))^2/slope,"(Subscript[π, N])"}},None},PlotRange->{{0,2(((maxp-cost)/(nfirms+1)))^2/slope},{0,100nfirms}}]
Show[SmoothHistogram[Flatten[simulatedquantities],{"Standardized", 1/4},PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{0,maxq/nfirms},Automatic},AxesOrigin->{0,0}]

Simulate Procedure 2 for N=16:

SeedRandom[1234]
Clear[nfirms,simulationburnin,simulationlength,simulatedquantities,simulatedtotals,simulatedprofits,meanprofits,dispprofits]
nfirms=16;
simulationburnin=nfirms^2;
simulationlength=20000; (* must be even *)
simulatedquantities=Module[{startcontinuous},
startcontinuous=With[{quantities=RandomReal[{0,maxq/nfirms},nfirms]},{quantities,1}];
Transpose[NestList[update2cour,startcontinuous,simulationburnin+simulationlength][[Range[simulationburnin+1,simulationburnin+simulationlength]]]]
][[1]];
simulatedtotals=Transpose@Table[Total[simulatedquantities,{2}],{nfirms}];
simulatedprofits=MapThread[(maxp-slope#2-cost)#1&,{simulatedquantities,simulatedtotals},2];
meanprofits=Map[Mean,simulatedprofits];
dispprofits=Map[RootMeanSquare,simulatedprofits-meanprofits];

Print["long run total = ",nfirms Mean[meanprofits]]
Print["benchmark total = ",N[nfirms(((maxp-cost)/(nfirms+1)))^2/slope]]
Print["long run volatility of the total = ",N[StandardDeviation[meanprofits]*nfirms]]
Print["correlation of mean and cross\[Hyphen]section SD = ",N[Correlation[meanprofits,dispprofits]]]
Print["correlation of mean and cross\[Hyphen]section COV = ",N[Correlation[meanprofits,dispprofits/meanprofits]]]
Print["long run COV = ",Mean[dispprofits/meanprofits]]

ListPlot[{Take[meanprofits,1000],Take[dispprofits,1000]},Joined->True,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{Automatic,{0}},AxesOrigin->{0,0},PlotRange->{0,2(((maxp-cost)/(nfirms+1)))^2/slope}]
SmoothHistogram[{meanprofits,dispprofits},{"Standardized", 1/4},Filling->Axis,PlotStyle->{{Darker@Blue},{Dashed,Darker@Red}},BaseStyle->Bold,ImageSize->250,Ticks->{{{(((maxp-cost)/(nfirms+1)))^2/slope,"(Subscript[π, N])"}},None},PlotRange->{{0,2(((maxp-cost)/(nfirms+1)))^2/slope},{0,100nfirms}}]
Show[SmoothHistogram[Flatten[simulatedquantities],{"Standardized", 1/4},PlotStyle->{Thick,Darker@Hue[.3,.5]}],BaseStyle->Bold,ImageSize->250,PlotRange->{{0,maxq/nfirms},Automatic},AxesOrigin->{0,0}]