11Value Functions and Entry Rules

We begin by implementing the computational algorithm to compute the equilibrium value functions that we present in the paper. The post-survival value function (StartMathJax)v_S(n,c)(StopMathJax) is computed recursively by iterating on a sequence of Bellman equations.

Recall the definitions of the entry thresholds in the paper,

The post-survival value function is given by

The above is the key equation that we will use to numerically compute the equilibrium. First, we consolidate the econometric error and obtain

Second, we invoke the distributional assumption on (StartMathJax)W(StopMathJax),

which gives us a closed form solution for the partial expectation,

where (StartMathJax)\Phi(\cdot)(StopMathJax) refers to the standard normal cumulative distribution function. The remaining two integrals in equation (3) can be expressed using the cumulative distribution function of (StartMathJax)W(StopMathJax):

We now implement the value function iteration on equation (3) in Matlab.

The function valueFunctionIteration requires the arguments Settings and Param. It returns the following four matrices:

• vS is a matrix of dimension (StartMathJax)(\check n + 1) \times \check c(StopMathJax) that stores the equilibrium post-survival value functions. By definition, the last row consists of zeros (and exists mostly for computational convenience).
• pEntry is a matrix of dimension (StartMathJax)(\check n + 1) \times \check c(StopMathJax) that stores the equilibrium entry probabilities, i.e. each element contains the probability that at least row firms are active post-entry under demand state column. Again, the last row consists of zeros. Thus, the (StartMathJax)n^{\text{th}}(StopMathJax) row of pEntry stores (StartMathJax)\Pr\left(w \lt \overline w_E(n, :)\right)(StopMathJax).
• pEntrySet is a matrix of dimension (StartMathJax)(\check n + 1) \times \check c(StopMathJax) that stores the equilibrium probabilities that exactly row firms are active post-entry under demand state column. Thus, the (StartMathJax)n^{\text{th}}(StopMathJax) row of pEntrySet stores (StartMathJax)\Pr\left(\overline w_E(n + 1, :) \leq w \lt \overline w_E(n, :)\right)(StopMathJax).
• pStay is a matrix of dimension (StartMathJax)\check n \times \check c(StopMathJax) that stores the equilibrium probabilities that row firms find certain survival profitable under demand state column. Thus, the (StartMathJax)n^{\text{th}}(StopMathJax) row of pStay stores (StartMathJax)\Pr\left(w \lt \overline w_S(n, :)\right)(StopMathJax).

The goal is to construct the (StartMathJax)(\check n + 1) \times \check c(StopMathJax) matrix vS, in which the last row is set to zero by construction. Each column in vS constitutes a fixed point to the Bellman equation characterized by equation (3). The procedure will follow a backward recursion on the number of firms, from (StartMathJax)\check n (StopMathJax) back to 1. For (StartMathJax)n = \check n(StopMathJax), the Bellman equation has (StartMathJax)v_S(\check n ,c)(StopMathJax) on both sides as entry cannot occur. For (StartMathJax)n \lt \check n(StopMathJax), the Bellman equation will depend only on (StartMathJax)v_S(n',c)(StopMathJax) for (StartMathJax)n' \gt n(StopMathJax), the values of which are already in memory because we are iterating backwards.

Allocate and intialize the post-survival value functions, the entry probabilities, and the probability of sure survival.

Now we begin the backward recursion, starting from nCheck. We will iterate on vS(n, :) using equation (3), which we map into the relevant Matlab variables below:

The expectation operator with respect to (StartMathJax)C'(StopMathJax) is implemented by left-multiplying the inside of the expectation in the equation above (a column vector with cCheck elements) by the transition probability matrix Param.demand.transMat. The iteration is complete when the difference vSdiff between vS(n, :) and its update n vSPrime does not exceed the stopping criterion, Settings.tolInner. Start by initializing vSdiff to 1 (which exceeds Settings.tolInner). We pre-compute (StartMathJax)\omega ^ 2(StopMathJax) and the demand grid (transposed) at the beginning, so we do not have to do so repeatedly inside the loops below.

Note that we only need to compute the entry probabilities outside of the value function iteration. We use the variable iter to keep track of the number of iterations for each value function iteration. Whenever iter exceeds Settings.maxIter, the value function iteration is terminated and a error is returned. This concludes valueFunctionIteration.

12Survival Rules

Our equilibrium computation algorithm valueFunctionIteration is fast because we do not need to completely characterize firms' survival strategies to compute the equilibrium value functions. We know that when firms mix between staying and exiting, the must receive a continuation value of zero. Firms use mixed strategies whenever the cost shock falls into the interval

which means that it is not profitable for all (StartMathJax)n(StopMathJax) active firms to continue, but the market is profitable enough for at least one firm to continue.

What we have not computed thus far is how firms mix between continuation and exit when the cost shock falls into this interval. The mixing probabilities (StartMathJax)a_S(StopMathJax) are implicitly defined by the indifference condition

The indifference condition states that when (StartMathJax)n(StopMathJax) active firms all use the survival rule a_S, the expected value from using survival rule a_S equals zero.

We compute the solution to (10) in mixingProbabilities.

This function computes the mixing probability for a market with N firms, demand state C, and cost shock W, where the post-survival value function is given by vS. The inputs N, C, and W are all of the same dimension.

The function returns the mixing probabilities aS which is of the same dimension as the inputs N, C, and W.

The function will solve (10) for (StartMathJax)a_S(StopMathJax) using Matlab's roots function. For the roots function to work, we need to transform (10) into polynomial form, which is given by

where the relevant Matlab variables are marked in bold font.

We allocate a vector of NaNs with the survival strategies that will subsequently be filled and then loop over each element in N.

Store the post-entry number of active firms in a scalar nE. Allocate the matrix matCoef to store the coefficients of the polynomial above. Then assemble the coefficients by looping from nE-1 to 0.

We then compute the candidate values for the mixing probabilities using roots, and nullify ([]) all values that are smaller than 0 or larger than 1. When the only root is really close to 0 or to 1, Matlab may return a couple of undesired complex roots as a bonus. We nullify these candidates as well. Finally, we pick the remaining root (which we know exists and is unique).

This concludes mixingProbabilities.

21Likelihood Step 1: Estimate (StartMathJax)\theta_C(StopMathJax)

This function computes the first step likelihood function from demand data. This function takes as inputs the Data structure (from which the matrix data.C will be used), the Settings structure, and the estimates vector which consists of (StartMathJax)\mu_C(StopMathJax) and (StartMathJax)\sigma_C(StopMathJax). The output of the function is the negative log likelihood ll and the (StartMathJax)(\check t -1)\cdot \check r \times1(StopMathJax) vector of likelihood contributions likeCont. Notice that the data contains (StartMathJax)\check t(StopMathJax) time periods, which gives us (StartMathJax)\check t -1(StopMathJax) transitions to construct the likelihood function.

We look for transitions from (StartMathJax)C_{t,r}(StopMathJax) to (StartMathJax)C_{t+1,r}(StopMathJax) for (StartMathJax)t=1,\ldots,\check t-1(StopMathJax). Start by constructing two matrices of size (StartMathJax)(\check t -1)\times \check r (StopMathJax) named from and to, which include (StartMathJax)C_{t,r}(StopMathJax) and (StartMathJax)C_{t+1,r}(StopMathJax), respectively, for (StartMathJax)t=1,\ldots,\check t-1(StopMathJax) and (StartMathJax)r=1,\ldots,\check r (StopMathJax). Thus, from and to are vectors with as many elements as there are transitions in the data.

We allocate the (StartMathJax)(\check t -1)\cdot \check r \times1(StopMathJax) likelihood contribution vector and set its value to NaN:

Assign muC and sigmaC, the values with respect to which we will maximize this likelihood function from the input estimates.

Now, compute the likelihood for each transition that is observed in the data given muC and sigmaC. For all transitions, calculate likelihood contributions according to Tauchen(1986). We make a distinction between transitions to interior points of the demand grid and transitions to points on the boundary.

Transitions to an interior point yield a likelihood contribution of

Similarly, transition to the lower and upper bound yield likelihood contributions of

and

respectively.

We then take the log of the contributions, sum them up, and return the negative log-likelihood contribution.

This concludes likelihoodStep1.

22Likelihood Step 2: Estimate (StartMathJax)(\theta_P,\theta_W)(StopMathJax)

We now construct the likelihood contributions that result from the number of firms evolving from (StartMathJax)n(StopMathJax) to (StartMathJax)n'(StopMathJax). Recall that we obtain cost-shock thresholds for entry and survival, defined by (StartMathJax)\overline w_{E}(n,c) \equiv \log v_{S}(n,c) - \log\left(1 + \varphi\right)(StopMathJax) and (StartMathJax)\overline w_{S}(n,c)\equiv \log v_{S}(n,c)(StopMathJax). We consider five mutually exclusive cases.

• Case I: (StartMathJax)\mathbf{n' \gt n}(StopMathJax). If the number of firms increases from (StartMathJax)n(StopMathJax) to (StartMathJax)n'(StopMathJax), then it must be profitable for the (StartMathJax)n'(StopMathJax)th firm to enter, but not for the (StartMathJax)(n'+1)(StopMathJax)th: (StartMathJax)\overline w_{E}(n'+1,c)\leq w \lt \overline w_{E}(n',c)(StopMathJax). The probability of this event is

  (23)

  (StartMathJax) \label{app:eq:llhcontr1} G_W\left[\overline w_{E}(n',c)\right]-G_W\left[\overline w_{E}(n'+1,c)\right]. (StopMathJax)

• Case II: (StartMathJax)\mathbf{0 \lt n' \lt n}(StopMathJax). If the number of firms decreases from (StartMathJax)n(StopMathJax) to (StartMathJax)n'(StopMathJax), with (StartMathJax)0 \lt n' \lt n(StopMathJax), then the realization of (StartMathJax)W(StopMathJax) must be such that firms exit with probability (StartMathJax)a_{S}(n,c,w)\in(0,1)(StopMathJax). Thus, the cost shock must be high enough so that (StartMathJax)n(StopMathJax) firms cannot survive profitably, (StartMathJax)w\geq\overline w_{S}(n,c)(StopMathJax), but low enough for a monopolist to survive profitably, (StartMathJax)w \lt \overline w_{S}(1,c)(StopMathJax). Given such value (StartMathJax)w(StopMathJax), (StartMathJax)N'(StopMathJax) is binomially distributed with success probability (StartMathJax)a_{S}(n,c,w)(StopMathJax) and population size (StartMathJax)n(StopMathJax). Hence, the probability of observing a transition from (StartMathJax)n(StopMathJax) to (StartMathJax)n'(StopMathJax) with (StartMathJax)0 \lt n' \lt n(StopMathJax) equals

  (24)

  (StartMathJax) \label{app:eq:llhcontr2} \int_{\overline w_{S}(n,c)}^{\overline w_{S}(1,c)} {n \choose n'} a_{S}(n,c,w)^{n'}\left[1-a_{S}(n,c,w)\right]^{n-n'} g_W(w)dw, (StopMathJax)

  where (StartMathJax)g_W(StopMathJax) is the density of (StartMathJax)G_W(StopMathJax). The integrand in (24) involves the mixing probabilities (StartMathJax)a_{S}(n,c,w)(StopMathJax). We discuss how we compute this integral in detail below.

• Case III: (StartMathJax)\mathbf{n'=0, n \gt 0}(StopMathJax). All firms exiting can be the result of two events. First, it is not profitable for even a single firm to continue, (StartMathJax)w \geq \overline w_{S}(1,c)(StopMathJax). Second, it is profitable for some but not all firms to continue, (StartMathJax)\overline w_{S}(n,c) \leq w \lt \overline w_{S}(1,c)(StopMathJax), firms exit with probability (StartMathJax)a_S(n,c)\in(0,1)(StopMathJax) as in Case II, and by chance none of the (StartMathJax)n(StopMathJax) firms survives. The probability of these events is

  (25)

  (StartMathJax) \label{app:eq:llhcontr3} 1-G_W\left[\overline w_{S}(1,c)\right] +\int_{\overline w_{S}(n,c)}^{\overline w_{S}(1,c)}\left[1-a_{S}(n,c,w)\right]^{n} g_W(w)dw. (StopMathJax)

• Case IV: (StartMathJax)\mathbf{n'=0, n=0}(StopMathJax). In this case, the market is populated by zero firms and it is not profitable for a monopolist to enter. The probability of this event is given by

  (26)

  (StartMathJax) \label{app:eq:llhcontr4} 1-G_W\left[\overline w_{E}(1,c)\right]. (StopMathJax)

• Case V: (StartMathJax)\mathbf{n' = n \gt 0}(StopMathJax). If there is neither entry nor exit, then either no firm finds it profitable to enter and all (StartMathJax)n(StopMathJax) incumbents find it profitable to stay, (StartMathJax)\overline w_{E}(n+1,c) \leq w \lt \overline w_{S}(n,c),(StopMathJax) or the (StartMathJax)n(StopMathJax) incumbents mix as in Cases II and III, but by chance end up all staying. The probability of these events is

  (27)

  (StartMathJax) \label{app:eq:llhcontr5} G_W\left[\overline w_{S}(n,c)\right]-G_W\left[\overline w_{E}(n+1,c)\right] + \int_{\overline w_{S}(n,c)}^{\overline w_{S}(1,c)} a_{S}(n,c,w)^{n} g_W(w)dw. (StopMathJax)

We compute the likelihood using the function likelihoodStep2 that requires as inputs the structures Data, Settings, and Param, and the vector estimates as inputs. It returns the scalar valued negative log-likelihood function ll and a column vector of length (StartMathJax)\check r \times (\check t - 1)(StopMathJax) containing the market-time-specific likelihood contributions, likeCont.

We start by mapping the vector estimates into the corresponding elements in the Param structure. We do this using anonymous functions that are defined in the structure Settings. By construction, Param.k is a vector of length (StartMathJax)\check{n}(StopMathJax). Param.phi and Param.omega are scalars.

Now we use valueFunctionIteration to solve the model by iterating on the post-survival value function. We also retrieve pStay, pEntry and pEntrySet, which are the probability of certain survival and the entry probabilities as described in valueFunctionIteration above.

Next we collect the transitions observed in the data and vectorize them. The column vectors from, to, and demand are all of length (StartMathJax)(\check t - 1) \times \check r(StopMathJax).

Here and throughout we will convert subscripts to linear indices using the Matlab function sub2ind.

We store the likelihood contributions in five vectors, each corresponding to one of the five cases outlined above. We allocate these vectors here and set all of their elements to zero.

Case I: We store all of the likelihood contributions resulting from entry in the vector llhContributionsCaseI.

Case II: We store all of the likelihood contributions resulting from exit to a non-zero number of firms in the vector llhContributionsCaseII.

Note that this case involves computing the integral over mixed strategy play, which we do in the function mixingIntegral. We document its content below.

Case III: We store all of the likelihood contributions resulting from transitions to zero (from a positive number of firms) in llhContributionsCaseIII.

Case IV: We store all of the likelihood contributions resulting from when the number of active firms remains at zero in llhContributionsCaseIV.

Case V: We store all of the likelihood contributions resulting from the number of firms staying the same in llhContributionsCaseV.

Finally, we sum up the likelihood contributions from the five cases and return the negative log likelihood function. When ll is not real valued, the negative log likelihood is set to inf.

If two outputs are requested, we also return the likelihood contributions:

This concludes likelihoodStep2.

We still need to specify what exactly happens in the function mixingIntegral.

Computing the integral resulting from mixed strategy play

For a given survival strategy (StartMathJax)a_S(n,c,w)(StopMathJax), the likelihood contribution from (purely) mixed strategy play is given by

The term inside the integral is the probability mass function of a binomial distribution function with success probability (StartMathJax)a_S(n,c,w)(StopMathJax). The survival strategies are defined by the indifference condition

In principle, we could compute the integral in (28) directly by naively numerically integrating over (StartMathJax)W(StopMathJax). In practice, it is computationally convenient to do a change of variables and integrate over the survival strategies (StartMathJax)a_S(n,c,\cdot)(StopMathJax) instead. To make an explicit distinction between the survival strategy (StartMathJax)a_S(n,c,w)(StopMathJax), which is a function of (StartMathJax)n(StopMathJax), (StartMathJax)c(StopMathJax), and (StartMathJax)w(StopMathJax), and the variable of integration, which is just a scalar. We will refer to the latter as (StartMathJax)p(StopMathJax). Thus, for a given value of (StartMathJax)p(StopMathJax), we need to find the value of (StartMathJax)w(StopMathJax) such that (StartMathJax)p = a_S(n,c,w)(StopMathJax).

Equation (29) defines the inverse (StartMathJax)a_S^{-1}(p;c,n)(StopMathJax) for which

This inverse function can be solved for analytically and it is given by

Then note that (StartMathJax)a_S^{-1}(1;c,n) = \log v_S(n,c)(StopMathJax) and (StartMathJax)a_S^{-1}(0;c,n) = \log v_S(1,c)(StopMathJax). We can write the likelihood contribution as an integral over (StartMathJax)p(StopMathJax):

where (StartMathJax)p_{1}, ..., p_{J}(StopMathJax) refer to the Gauss-Legendre nodes and (StartMathJax)w_{1}, ..., w_{J}(StopMathJax) to the corresponding weights. Notice that the integration bounds are now 0 and 1 since if (StartMathJax)w \lt \log v_S(n,c)(StopMathJax) the firms surely survive and when (StartMathJax)w \gt \log v_S(1,c)(StopMathJax) the firms surely exit. Differentiation of (StartMathJax)a_S^{-1}(p;c,n)(StopMathJax) gives

Now, compute the matrix mixingDensity using (34). mixingDensity is of dimension Settings.integrationLength by Settings.cCheck by Settings.nCheck. It is defined as

The element (StartMathJax)(p_{j}, c, n)(StopMathJax) gives us the density of the mixing probability (StartMathJax)p_{j}(StopMathJax) when demand equals (StartMathJax)c(StopMathJax) and the current number of incumbents is (StartMathJax)n(StopMathJax).

In the function mixingIntegral we compute the integral in equation (28) for a range of different combinations of (StartMathJax)n(StopMathJax), (StartMathJax)n'(StopMathJax), and (StartMathJax)c(StopMathJax) using the change of variable introduced above. The function mixingIntegral takes as arguments the vectors from, to, and demand, which correspond to (StartMathJax)n(StopMathJax), (StartMathJax)n'(StopMathJax), and (StartMathJax)c(StopMathJax), respectively. The function also takes as argument vS, the equilibrium post-survival value functions, and the Param and Settings structures. The function returns a vector llhContributionsMixing that is of the same dimension as the inputs from, to, and demand.

Note that mixed strategy play is only relevant for markets with at least two firms. We define the auxiliary variable expaSInv, which equals expaSInv = exp(aSinv(:, :, n)). dexpaSInvdP is another auxiliary variable that stores the derivative of expaSInv with respect to p. Then assemble expaSInv and dexpaSInvdP by looping over the possible outcomes from mixing nPrime=0,...,n. nPrime=0 refers to the case when all firms leave due to mixed strategy play and nPrime=n refers to the case when all firm stay due to mixed strategy play. Then, use val to compute aSinv and compute the derivative with respect to p and store it as daSinvdP. Lastly, we compute the mixing density.

We pre-compute nchoosekMatrixPlusOne which is a matrix of size (StartMathJax)\check n + 1(StopMathJax) by (StartMathJax)\check n + 1(StopMathJax), where element (StartMathJax)(i,j)(StopMathJax) contains (StartMathJax)i - 1 \choose j - 1(StopMathJax). The copious naming and indexing convention is owed to the fact that Matlab indexing starts at one, not zero, so element (StartMathJax)(1,1)(StopMathJax) corresponds to (StartMathJax)0 \choose 0(StopMathJax). Pre-computing this matrix is helpful, because factorial operations are computationally demanding.

We then set up to compute the matrix mixingDensity, as given by (34).

With the matrix mixingDensity in hand, we can compute the likelihood contributions from mixed strategy play using (32). Note that this time, we cannot avoid the use of loops altogether. However, we only need to loop over all those observations where “purely” mixed strategy play does in fact occur.

Note that it is possible to improve the performance of the code by only calling this function once instead of three times as we currently do in the likelihood function computation. However, that makes the code somewhat more difficult to follow, which is why we opted for the slightly slower version.

23Likelihood Step 3: Estimate (StartMathJax)(\theta_C, \theta_P, \theta_W)(StopMathJax)

This function computes the full information likelihood. This step only involves taking the product of the likelihood contributions from the first two steps. Here we will also go over the computation of standard errors. The function likelihoodStep3 requires the structures Data, Settings, and Param, and the vector estimates as inputs. The arguments returned include the negative log likelihood function ll, a vector of standard errors se, a vector of the likelihood contributions likeCont, and the covariance matrix covariance.

Create the transition probability matrix based on the estimates for muC and sigmaC which are stored as the last two elements in estimates. Then retrieve the vectors of likelihood contributions from the first two steps and compute the negative log likelihood function of the third step:

This completes the construction of the full information negative log likelihood function. When only one output argument is requested, then the function is done.

When three output arguments are requested, the function returns the likelihood contributions, which are simply the element-by-element product of the likelihood contributions vectors from the first two steps. In this case, we return no standard errors, i.e. we set se=[].

Now, consider the case, when exactly two output arguments are requested. In this case, we want to compute the standard errors. As discussed in the paper, standard errors are computed using the outer-product-of-the-gradient estimator of the information matrix. When two output arguments are requested when calling likelihootStep3.m, the function will return standard errors in addition to the log likelihood function. When three output arguments are requested, the function also returns the likelihood contributions. Define the matrix of perturbations, which is simply a diagonal matrix with Settings.fdStep on each diagonal element.

Next, get the likelihood contributions at the parameter values in estimates:

Now, given the likelihood contribution (StartMathJax)\ell(\theta) \equiv \ell(\theta_j, \theta_{-j})(StopMathJax) we compute for each parameter (StartMathJax)\theta_j(StopMathJax) the positively and negatively perturbed likelihood contributions (StartMathJax)\ell(\theta_j+\epsilon,\theta_{-j})(StopMathJax) and (StartMathJax)\ell(\theta_j-\epsilon,\theta_{-j})(StopMathJax). The gradients of the negative log likelihood contributions are then computed using central finite differences:

The matrix of gradient contributions gradCont has (StartMathJax)(\check t - 1)\cdot \check r (StopMathJax) rows and one column for each parameter with respect to which we are differentiating the logged likelihood contributions.

We now have the matrix gradCont where each column is the score with respect to an estimable parameter. This matrix is used to compute the Hessian (full information matrix). We take the sum of the outer-product of the market specific gradients over all markets. A market specific gradient is a row in gradCont. Looping over all rows, in each iteration we compute the outer product of a market specific gradient and add them all up:

The covariance matrix can be obtained from inverting the Hessian. The standard errors are the square roots of the diagonal entries of the covariance matrix.

This concludes likelihoodStep3.

31Draw Number of Firms

The function randomFirms randomly draws a realization of the surviving number of active firms (StartMathJax)N'(StopMathJax) given the current period's initial number of active firms (StartMathJax)n(StopMathJax), demand state (StartMathJax)c(StopMathJax) and cost shock (StartMathJax)w(StopMathJax). The computation of the number of firms follows the construction of the unique symmetric Nash equilibrium of the one-shot survival game, which takes advantage of the monotonicity of (StartMathJax)v_S(n,c)(StopMathJax). Given (StartMathJax)n(StopMathJax), (StartMathJax)c(StopMathJax), and (StartMathJax)w(StopMathJax), the possible realizations of (StartMathJax)N'(StopMathJax) can be classified using the following cases:

• (StartMathJax)v_S(1,c) \leq \exp (w)(StopMathJax): In this case, even a monopolist would not be willing to remain active. Thus, all firms leave the market for sure and there is no entry. This case is only relevant when (StartMathJax)n \gt 0(StopMathJax), i.e. when there is a positive number of incumbents. If the number of incumbents is equal to zero, and the above condition holds, then there will not be any entry and the number of firms remains at zero.
• (StartMathJax)v_S(n,c) \leq \exp (w) \lt v_S(1,c)(StopMathJax): In this case, a monopolist finds survival profitable, but the current number of firms (StartMathJax)n(StopMathJax) does not. Thus, some firms will be leaving the market and firms' survival strategies satisfy (StartMathJax)a_S(n,c, w) \in (0,1)(StopMathJax). (StartMathJax)N'(StopMathJax) follows a binomial distribution with (StartMathJax)n(StopMathJax) trials and success probability (StartMathJax)a_S(n,c,w)(StopMathJax). This case is only well defined when (StartMathJax)n \gt 1(StopMathJax).
• (StartMathJax)v_S(n,c) \gt \exp (w) (StopMathJax): All (StartMathJax)n(StopMathJax) incumbents find survival profitable. In addition, there may be some entry. We will use the entry strategies to compute (StartMathJax)n'(StopMathJax). In fact, we will count for how many (StartMathJax)n'(StopMathJax) the entry condition (StartMathJax)v_S(n',c) \gt (1+\varphi) \exp (w) (StopMathJax) is satisfied. This case is only relevant if (StartMathJax)n \lt \check n(StopMathJax), i.e. the current number of incumbents is strictly less than (StartMathJax)\check n(StopMathJax). If (StartMathJax)n=\check n(StopMathJax), there will not be any entry by construction.

The function takes as inputs last period's post-survival number of active firms N, the current number of consumers C, the realized cost shock W, and the equilibrium post-survival value function, vS. N, C, and W are column vectors of length rCheck containing one entry per market. The output is the column vector Nprime of length rCheck, which contains the post-survival number of active firms in each market.

This concludes randomFirms.

32Assemble Data Set

This function generates data (StartMathJax)(N,C,W)(StopMathJax) for Settings.rCheck markets and Settings.tCheck time periods. Demand in the first period is drawn from the ergodic distribution. To eliminate initial condition effects, Settings.tBurn+Settings.tCheck periods will be generated, but only the last Settings.tCheck periods will be used. The function requires the structures Settings and Param as inputs. The function returns a Data structure with three (StartMathJax)\check t \times \check r(StopMathJax) matrices N, C, and W, where each entry corresponds to some time period (StartMathJax)t(StopMathJax) in some market (StartMathJax)r(StopMathJax). N contains the number of active firms. C contains the index for the demand grid logGrid that corresponds to the demand state. W consists of realizations of the cost shocks. In practice, W is unobserved to the econometrician and thus we will not use it during the estimation. Here, we return W as additional information that can be used to debug the program or to compute additional statistics, such as the average realized fixed costs or the average realized entry costs.

We now compute the post-survival equilibrium value functions using the function valueFunctionIteration:

The cost shocks are iid across markets and periods, so we can draw them all at once and store them in the matrix W which is of dimension Settings.tBurn + Settings.tCheck by Settings.rCheck.

Next, we draw an initial demand state from the ergodic distribution of the demand process for each market using the randomDiscr function, which is documented in the appendix. With the initial demand state for each market in hand, for each time period we use the transition matrix to draw the current demand state given the previous state:

The initial number of firms in each market is drawn randomly from a discrete uniform distribution on (StartMathJax)\{{1,2,\ldots ,\check n }\}(StopMathJax). In each period following the first one, the number of firms is generated using the randomFirms function, which randomly draws realizations of the cost shocks and then uses the firms' equilibrium strategies to update the number of active firms. We draw the numbers of firms separately for the burn-in phase and the real sample, because for the latter we also want to store the realized fixed and entry costs.

Now we have the matrices N, C, and W; each of dimension Settings.tBurn + Settings.tCheck by Settings.rCheck. From this, we store only the last Settings.tCheck periods in the structure Data:

This concludes the function dgp.

61Compute Markov Process

This function computes the transition probability matrix transMat and the ergodic distribution ergDist of the demand process. If this function is called using only Param and Settings as inputs, then Param.demand.transMat is computed using the true values of the demand process parameters (StartMathJax)(\mu_C, \sigma_C)(StopMathJax) that are stored in Param.truth.step1. This is used in example.m when we generate synthetic data on the number of consumers in each market. In this case, we will also compute the ergodic distribution Param.demand.ergDist. Alternatively, if the function is called with the additional inputs muC and sigmaC, then these parameter values will be used to compute Param.demand.transMat. In this case, Param.demand.ergDist will not be computed.

The code starts with extracting the relevant variables from the Settings structure for convenience. If only two input arguments are passed to the function, then the true values for (StartMathJax)(\mu_C, \sigma_C)(StopMathJax) are used.

The transition probability matrix is computed according to the method outlined by Tauchen(1986), which is used to discretize a continuous (and bounded) stochastic process. We use the standard convention for transition probability matrices. That is, for a transition probability matrix (StartMathJax)\Pi(StopMathJax), each entry (StartMathJax)\Pi_{i,j}(StopMathJax) gives the probability of transitioning from state (StartMathJax)i(StopMathJax) (row) to state (StartMathJax)j(StopMathJax) (column). The transition matrix is of dimension (StartMathJax)\check c \times \check c(StopMathJax). The idea of the Tauchen method is intuitive - we assumed the growth of (StartMathJax)C(StopMathJax) to be normally distributed with parameters muC and sigmaC. Since this is a continuous distribution, while the state space is discrete, we treat a transition to a neighborhood around (StartMathJax)c_{[j]}(StopMathJax) as a transition to (StartMathJax)c_{[j]}(StopMathJax) itself. These neighborhoods span from one mid-point between two nodes on logGrid to the next mid-point, i.e. (StartMathJax)\left[\log c_{[j]}-\frac{d}{2},\log c_{[j]}+\frac{d}{2}\right](StopMathJax). Transitions to the end-points of logGrid follow the same logic. Distinguishing between transitions to interior points ((StartMathJax)j=2,\ldots,\check c -1(StopMathJax)) and transitions to end-points ((StartMathJax)j=1(StopMathJax) or (StartMathJax)j=\check c(StopMathJax)), we have three cases.

This completes the construction of the transition matrix transMat, which we now store in the Param structure.

Next, we compute the ergodic distribution of the demand process ergDist. The ergodic distribution is only computed for the true parameter values, i.e. when the number of input arguments equals two. The ergodic distribution is the eigenvector of the transpose of the transition matrix with eigenvalue equal to 1 after normalizing it such that its entries sum to unity. The Perron-Frobenius Theorem guarantees that such an eigenvector exists and that all eigenvalues are not greater than 1 in absolute value. We store the eigenvectors of the transition matrix as columns in eigenVecs, in decreasing order of eigenvalues from left to right. The eigenvector with unit eigenvalue is normalized by dividing each element by the sum of elements in the vector, so that it sums to 1 and thus is the ergodic distribution, stored as ergDist. Finally, we store ergDist in the Param structure.

We conclude by checking for two numerical problems that may arise. Firstly, confirm that the greatest eigenvalue is sufficiently close to 1 and return an error if it is not. Secondly, due to approximation errors, it may be the case that not all elements in the eigenvector are of the same sign (one of them may be just under or above zero). This will undesirably result in an ergodic distribution with a negative entry.

62Draw from Discrete Distribution

This function randomly draws realizations from discrete probability distributions using the inverse cumulative distribution function method. The idea is to first draw a realization from a uniform distribution and then map that draw, a probability, into the corresponding mass point of the discrete distribution. This function is written to draw one realization each from K distributions at a time. Each of the K distributions is assumed to be defined over M points.

The function takes as input a matrix P which is of dimension M times K. In this matrix, each column corresponds to a probability mass function.

The result is a vector iX of length K with realizations from the discrete distributions described by P.

63Compute Gauss-Legendre Weights

This function is for computing definite integrals using Legendre-Gauss Quadrature. It computes the Legendre-Gauss nodes and weights on an interval [a,b] with number of nodes N.

Suppose you have a continuous function (StartMathJax)f(x)(StopMathJax) which is defined on [a,b] which you can evaluate at any (StartMathJax)x(StopMathJax) in [a,b]. Simply evaluate it at all of the values contained in the x vector to obtain a vector (StartMathJax)f(x)(StopMathJax). Then compute the definite integral using sum(f .* w);

• Abbring, J. H., J. R. Campbell, J. Tilly, N. Yang (2018): "Very Simple Markov-Perfect Industry Dynamics: Theory," Econometrica, 86, 721-735.
• Abbring, J. H., J. R. Campbell, J. Tilly, N. Yang (2017): "Very Simple Markov-Perfect Industry Dynamics: Empirics," mimeo.
• Dube, J.-P., J. T. Fox, and C.-L. Su (2012): "Improving the Numerical Performance of Static and Dynamic Aggregate Discrete Choice Random Coefficients Demand Estimation," Econometrica, 80, 2231-2267.
• Rust, J. (1987): "Optimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcher," Econometrica, 55, 999-1033.
• Tauchen, G. (1986): "Finite State Markov-Chain Approximations to Univariate and Vector Autoregressions," Economic Letters, 20, 177-181.