11Flow Payoffs

The function flowpayoffs computes the mean (over (StartMathJax)\varepsilon(StopMathJax)) flow payoffs, (StartMathJax)u_0(x,a)(StopMathJax) and (StartMathJax)u_1(x,a)(StopMathJax), for each profit and past choice pair (StartMathJax)(x,a)\in{\cal X}\times\{0,1\}(StopMathJax).

It requires the following input arguments:

{supportX} a (StartMathJax)K\times 1(StopMathJax) vector with the support points of the profit state (StartMathJax)X_t(StopMathJax) (the elements of (StartMathJax)\cal{X}(StopMathJax), consistently ordered with the Markov transition matrix (StartMathJax)\Pi(StopMathJax));

{beta} a (StartMathJax)2\times 1(StopMathJax) vector that contains the intercept ((StartMathJax)\beta_0(StopMathJax)) and profit state slope ((StartMathJax)\beta_1(StopMathJax)) of the net payoffs to choice (StartMathJax)1(StopMathJax); and

{delta} a (StartMathJax)2\times 1(StopMathJax) vector that contains the firm's exit ((StartMathJax)\delta_0(StopMathJax)) and entry ((StartMathJax)\delta_1(StopMathJax)) costs.

It returns

{u0} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)u_0(x^i,j-1)(StopMathJax) and

{u1} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)u_1(x^i,j-1)(StopMathJax) .

That is, the rows correspond to the support points of (StartMathJax)X_t(StopMathJax), and the columns to the choice in the previous period, (StartMathJax)A_{t-1}(StopMathJax).

The function flowpayoffs first stores the number (StartMathJax)K(StopMathJax) of elements of supportX in a scalar nSuppX.

Then, it constructs a (StartMathJax)K\times 2(StopMathJax) matrix u0 with the value of

and a (StartMathJax)K\times 2(StopMathJax) matrix u1 with the value of

You can change the specification of the flow profits by editing these two lines of code.

12Bellman Operator

The expected discounted profits net of (StartMathJax)\varepsilon_{t}(a)(StopMathJax) immediately following choice (StartMathJax)a(StopMathJax) at time (StartMathJax)t(StopMathJax), (StartMathJax)U_a(X_t,A_{t-1})(StopMathJax), satisfy a recursive system (StartMathJax)U_0=\Psi_0(U_0,U_1)(StopMathJax) and (StartMathJax)U_1=\Psi_1(U_0,U_1)(StopMathJax) or, with (StartMathJax)U\equiv(U_0,U_1)(StopMathJax) and (StartMathJax)\Psi\equiv(\Psi_0,\Psi_1)(StopMathJax), simply (StartMathJax)U=\Psi(U)(StopMathJax) (see e.g. Rust (1994)). Here, (StartMathJax)\Psi(StopMathJax) is a Bellman-like operator that embodies the model's dynamic specification and that depends on the flow payoffs (StartMathJax)u_0(StopMathJax) and (StartMathJax)u_1(StopMathJax), the Markov transition matrix (StartMathJax)\Pi(StopMathJax), and the discount factor (StartMathJax)\rho(StopMathJax). Its elements are the operators (StartMathJax)\Psi_0(StopMathJax) and (StartMathJax)\Psi_1(StopMathJax), with (StartMathJax)\Psi_a(StopMathJax) implicitly defined by the right hand side of

where (StartMathJax)R_a(x)\equiv\mathbb{E}\left[\max\{U_0(x,a)+\varepsilon_{t+1}(0),U_1(x,a)+\varepsilon_{t+1}(1)\}\right](StopMathJax) is McFadden (1981)'s social surplus for the binary choice problem with utilities (StartMathJax)U_0(x,a)+\varepsilon_{t+1}(0)(StopMathJax) and (StartMathJax)U_1(x,a)+\varepsilon_{t+1}(1)(StopMathJax). The first term in the right-hand side of (1) equals period (StartMathJax)t(StopMathJax)'s flow profits following choice (StartMathJax)a(StopMathJax), net of (StartMathJax)\varepsilon_t(a)(StopMathJax). The second term equals the expected discounted profits from continuing into period (StartMathJax)t+1(StopMathJax) with choice (StartMathJax)a(StopMathJax) (note that, given the choice (StartMathJax)a(StopMathJax), these continuation payoffs do not depend on the past choice (StartMathJax)A_{t-1}(StopMathJax)). The (mean-zero) extreme value assumptions on (StartMathJax)\varepsilon_{t+1}(0)(StopMathJax) and (StartMathJax)\varepsilon_{t+1}(1)(StopMathJax) imply that

The function bellman applies the operator (StartMathJax)\Psi(StopMathJax) once to input values of (StartMathJax)U(StopMathJax) and returns (StartMathJax)\Psi(U)(StopMathJax).

It requires the following input arguments:

{capU0} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)U_0(x^i,j-1)(StopMathJax);

{capU1} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)U_1(x^i,j-1)(StopMathJax);

{u0} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)u_0(x^i,j-1)(StopMathJax);

{u1} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)u_1(x^i,j-1)(StopMathJax);

{capPi} the (StartMathJax)K\times K(StopMathJax) Markov transition matrix (StartMathJax)\Pi(StopMathJax) for (StartMathJax)\{X_t\}(StopMathJax), with typical element (StartMathJax)\Pi_{ij}=\Pr(X_{t+1}=x^j|X_t=x^i)(StopMathJax); and

{rho} a scalar with the value of the discount factor (StartMathJax)\rho(StopMathJax).

It returns

{capU0} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)U_0(x^i,j-1)(StopMathJax);

{capU1} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)U_1(x^i,j-1)(StopMathJax).

To this end, bellman first computes the surpluses (StartMathJax)R_0(x)(StopMathJax) and (StartMathJax)R_1(x)(StopMathJax) in (2) for all (StartMathJax)x\in\cal{X}(StopMathJax) and stacks these in (StartMathJax)K\times 1(StopMathJax) vectors r0 and r1.

Then, it applies (1) to compute new values of capU0 and capU1.

Here, the conditional expectation over (StartMathJax)X_{t+1}(StopMathJax) in (1) is taken by premultiplying the vectors r0 and r1 by the Markov transition matrix capPi. The vectors r0 and r1 are postmultiplied by [1 1] because the surpluses, and therefore the continuation payoffs, are independent of the past choice that indexes the columns of capU0 and capU1.

The logit assumption only affects the operator (StartMathJax)\Psi(StopMathJax), and therefore the function bellman, through the specification of the surpluses (StartMathJax)R_0(x)(StopMathJax) and (StartMathJax)R_1(x)(StopMathJax) in (2). If you want to change the logit assumption, you should change the computation of r0 and r1 (and make sure to adapt the computation of choice probabilities and inverse choice probabilities elsewhere as well).

41First Stage: State Transitions

The function estimatePi computes a frequency estimate of the Markov transition matrix (StartMathJax)\Pi(StopMathJax) from state transition data.

It requires the following input arguments:

{iX} a (StartMathJax)T\times N(StopMathJax) matrix with indices of observed states (StartMathJax)x_{tn}(StopMathJax) in (StartMathJax){\cal X}(StopMathJax) (for example, if (StartMathJax)x_{11}=x^3(StopMathJax), then the first element of iX is 3, not (StartMathJax)x^3(StopMathJax));

{nSuppX} the scalar number (StartMathJax)K(StopMathJax) of support points of the profit state (StartMathJax)X_t(StopMathJax) (the number of elements of (StartMathJax)\cal{X}(StopMathJax))

and returns

{piHat} an estimate (StartMathJax)\hat\Pi(StopMathJax) of the (StartMathJax)K\times K(StopMathJax) Markov transition matrix (StartMathJax)\Pi(StopMathJax) for (StartMathJax)\{X_t\}(StopMathJax), with typical element (StartMathJax)\hat\Pi_{ij}(StopMathJax) equal to the sample frequency of transitions to state (StartMathJax)j(StopMathJax) among all transitions from state (StartMathJax)i(StopMathJax).

The function estimatePi first stores the number of time periods (StartMathJax)T(StopMathJax) in a scalar nPeriods.

Then, for each pair (StartMathJax)(i,j)\in\{1,\ldots,K\}\times\{1,\ldots,K\}(StopMathJax), it estimates the probability (StartMathJax)\Pi_{ij}=\Pr(X_{t+1}=x^j|X_t=x^i)(StopMathJax) by the appropriate sample frequency, the number of transitions from (StartMathJax)i(StopMathJax) to (StartMathJax)j(StopMathJax) divided by the total number of transitions from (StartMathJax)i(StopMathJax) in the data iX.

Note that estimatePi requires a positive number of transition observations from each state. More generally, the frequency estimator that it implements only performs well with samples that are large relative to the state space. With relatively small samples, the frequency estimator should be replaced by one that smoothes across support points.

42Second Stage: Choices

The function negLogLik computes minus the log partial likelihood for the conditional choice part of the data. Optionally, it also returns minus the corresponding score vector and an estimate of the information matrix for the parameter (sub)vector (StartMathJax)\theta\equiv(\beta_0,\beta_1,\delta_1)'(StopMathJax) (the scores are specific to the estimation example in Section 5's script and should be adapted for inference on other parameters).

The function negLogLik requires the following input arguments:

{choices} a (StartMathJax)T\times N(StopMathJax) matrix with choice observations (StartMathJax)a_{tn}(StopMathJax);

{iX} a (StartMathJax)T\times N(StopMathJax) matrix with indices of observed states (StartMathJax)x_{tn}(StopMathJax) in (StartMathJax){\cal X}(StopMathJax) (for example, if (StartMathJax)x_{11}=x^3(StopMathJax), then the first element of iX is 3, not (StartMathJax)x^3(StopMathJax));

{supportX} a (StartMathJax)K\times 1(StopMathJax) vector with the support points of the profit state (StartMathJax)X_t(StopMathJax) (the elements of (StartMathJax)\cal{X}(StopMathJax), consistently ordered with the Markov transition matrix (StartMathJax)\Pi(StopMathJax));

{capPi} the (possibly estimated) (StartMathJax)K\times K(StopMathJax) Markov transition matrix (StartMathJax)\Pi(StopMathJax) for (StartMathJax)\{X_t\}(StopMathJax), with typical element (StartMathJax)\Pi_{ij}=\Pr(X_{t+1}=x^j|X_t=x^i)(StopMathJax);

{beta} a (StartMathJax)2\times 1(StopMathJax) vector that contains the intercept ((StartMathJax)\beta_0(StopMathJax)) and profit state slope ((StartMathJax)\beta_1(StopMathJax)) of the flow payoffs to choice (StartMathJax)1(StopMathJax);

{delta} a (StartMathJax)2\times 1(StopMathJax) vector that contains the firm's exit ((StartMathJax)\delta_0(StopMathJax)) and entry ((StartMathJax)\delta_1(StopMathJax)) costs;

{rho} a scalar with the value of the discount factor (StartMathJax)\rho(StopMathJax);

{flowpayoffs} a handle of a function [u0,u1]=flowpayoffs(supportX,beta,delta) that computes the mean flow payoffs (StartMathJax)u_0(StopMathJax) and (StartMathJax)u_1(StopMathJax);

{bellman} a handle of a function [capU0,capU1] = bellman(capU0,capU1,u0,u1,capPi,rho) that iterates once on (StartMathJax)\Psi(StopMathJax);

{fixedPoint} a handle of a function [capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,bellman,capU0,capU1) that computes the fixed point (StartMathJax)U(StopMathJax) of (StartMathJax)\Psi(StopMathJax); and

{tolFixedPoint} a scalar tolerance level that is used to determine convergence of the successive approximations of the fixed point (StartMathJax)U(StopMathJax) of (StartMathJax)\Psi(StopMathJax).

It returns

{nll} a scalar with minus the log partial likelihood for the conditional choices

and optionally

{negScore} a (StartMathJax)3\times 1(StopMathJax) vector with minus the partial likelihood score for (StartMathJax)\theta(StopMathJax) and

{informationMatrix} a (StartMathJax)3\times 3(StopMathJax) matrix with the sum of the (StartMathJax)N(StopMathJax) outer products of the individual contributions to the score for (StartMathJax)\theta(StopMathJax).

The function negLogLik first stores the number (StartMathJax)K(StopMathJax) of elements of supportX in a scalar nSuppX.

Next, it computes the flow payoffs (StartMathJax)u_0(StopMathJax) (u0) and (StartMathJax)u_1(StopMathJax) (u1), the choice-specific net expected discounted values (StartMathJax)U_0(StopMathJax) (capU0) and (StartMathJax)U_1(StopMathJax) (capU1), their contrast (StartMathJax)\Delta U(StopMathJax) (deltaU), and the implied probabilities (StartMathJax)1/\left[1+\exp(\Delta U)\right](StopMathJax) of not serving the market (pExit) for the inputted parameter values. Note that this implements the NFXP procedure's inner loop.

Log Partial Likelihood

The contribution to the likelihood of firm (StartMathJax)n(StopMathJax)'s choice in period (StartMathJax)t(StopMathJax) is the conditional choice probability

with (StartMathJax)a_{0n}=0(StopMathJax). The function negLogLik first computes these probabilities for each firm (StartMathJax)n(StopMathJax) and period (StartMathJax)t(StopMathJax) and stores them in a (StartMathJax)T\times N(StopMathJax) matrix p. Then, it returns minus the sum of their logs, the log partial likelihood for the conditional choices, in nll.

Score

If two or more output arguments are demanded from negLogLik, it computes and returns minus the partial likelihood score for (StartMathJax)\theta(StopMathJax) (the derivative of minus the log partial likelihood with respect to (StartMathJax)\theta(StopMathJax)), in negScore.

Firm (StartMathJax)n(StopMathJax)'s contribution to the score equals

Its calculation requires that we compute (StartMathJax)\partial\Delta U/\partial\theta(StopMathJax). Recall that (StartMathJax)U(StopMathJax), and therewith (StartMathJax)\Delta U(StopMathJax), is only implicitly given by (StartMathJax)U=\Psi(U)(StopMathJax). Note that (StartMathJax)U=(U_0,U_1)(StopMathJax) is defined on a set with (StartMathJax)2K(StopMathJax) points, so that (StartMathJax)U(StopMathJax) can be represented by a (StartMathJax)4K\times 1(StopMathJax) vector and (StartMathJax)\Psi(StopMathJax) by a mapping from (StartMathJax)\mathbb{R}^{4K}(StopMathJax) into (StartMathJax)\mathbb{R}^{4K}(StopMathJax). Specifically, (StartMathJax)U(StopMathJax) can be represented by the (StartMathJax)4K\times 1(StopMathJax) vector (StartMathJax)\bar U(StopMathJax) that lists the values of (StartMathJax)U_0(StopMathJax) and (StartMathJax)U_1(StopMathJax) on their domain,

and that satisfies

with (StartMathJax)\tilde\Psi_\theta:\mathbb{R}^{4K}\rightarrow\mathbb{R}^{4K}(StopMathJax) an appropriately rearranged version of (StartMathJax)\Psi(StopMathJax) (note that we made its dependence on (StartMathJax)\theta(StopMathJax) explicit). With this alternative representation of (StartMathJax)U(StopMathJax) and (StartMathJax)\Psi(StopMathJax) in place, we can solve

for (StartMathJax)\partial\bar U/\partial\theta'(StopMathJax) (see also Rust (1994), p.3110), where (StartMathJax)I_i(StopMathJax) is a (StartMathJax)i\times i(StopMathJax) unit matrix;

with

and (StartMathJax)O_i(StopMathJax) is a (StartMathJax)i\times i(StopMathJax) matrix of zeros; and

The function negLogLik first computes (StartMathJax)D_{00}(StopMathJax) (d00), (StartMathJax)D_{01}(StopMathJax) (d01), (StartMathJax)D_{10}(StopMathJax) (d10), (StartMathJax)D_{11}(StopMathJax) (d11) and constructs (StartMathJax)\partial \tilde\Psi_\theta(\bar U)/\partial\bar U'(StopMathJax) (dPsi_dUbar).

It then computes (StartMathJax)\partial\tilde\Psi_\theta(\bar U)/\partial\theta'(StopMathJax) (dPsi_dTheta; note that the following line is the only code in negLogLik that needs to be changed if other parameters than those in (StartMathJax)\theta(StopMathJax) are estimated).

Next, it computes (StartMathJax)\partial\bar U/\partial\theta'(StopMathJax) (dUbar_dTheta) and (StartMathJax)\partial\Delta U/\partial\theta'(StopMathJax) (dDeltaU_dTheta).

Finally, it computes the (StartMathJax)1\times 3(StopMathJax) vector (StartMathJax)-\partial\log\left[\prod_{t=1}^T p(a_{tn}|x_{tn},a_{(t-1)n})\right]/\partial \theta'(StopMathJax) for each (StartMathJax)n(StopMathJax), stacks these individual (minus) score contributions in the (StartMathJax)N\times 3(StopMathJax) matrix negFirmScores, and sums them to compute minus the score vector, negScore.

Information Matrix

If three output arguments are demanded, negLogLik computes the sum of the (StartMathJax)N(StopMathJax) outer products of the individual score contributions,

and returns it in informationMatrix.

When evaluated at an estimate of (StartMathJax)\theta(StopMathJax), this provides an estimate of the expected partial likelihood information matrix for (StartMathJax)\theta(StopMathJax). This estimate can in turn be used to estimate the variance-covariance matrix, and thus the standard errors, of the maximum partial likelihood estimator (StartMathJax)\hat\theta(StopMathJax) of (StartMathJax)\theta(StopMathJax).

51Simulating Data

First, we set the number of time periods (nPeriods) and firms (nFirms) that we would like to have in our sample.

We also set the tolerance tolFixedPoint on the fixed point (StartMathJax)U(StopMathJax) of (StartMathJax)\Psi(StopMathJax) that we will use to determine the simulation's entry and exit rules. This same tolerance will also be used when solving the model in the inner loop of the NFXP procedure.

Next, we specify the values of the model's parameters used in the simulation:

{nSuppX} the scalar number (StartMathJax)K(StopMathJax) of elements of (StartMathJax){\cal X}(StopMathJax);

{supportX} the (StartMathJax)K\times 1(StopMathJax) vector (StartMathJax){\cal X}(StopMathJax) with the support points of (StartMathJax)X_t(StopMathJax);

{capPi} the (StartMathJax)K\times K(StopMathJax) Markov transition matrix (StartMathJax)\Pi(StopMathJax) for (StartMathJax)\{X_t\}(StopMathJax), with typical element (StartMathJax)\Pi_{ij}=\Pr(X_{t+1}=x^j|X_t=x^i)(StopMathJax);

{beta} the (StartMathJax)2\times 1(StopMathJax) vector (StartMathJax)\beta(StopMathJax) with the parameters of the flow profit of active firms;

{delta} the (StartMathJax)2\times 1(StopMathJax) vector of exit and entry costs (StartMathJax)\delta(StopMathJax); and

{rho} the scalar discount factor (StartMathJax)\rho(StopMathJax).

For these parameter values, we compute the flow payoffs (StartMathJax)u_0(StopMathJax) (u0) and (StartMathJax)u_1(StopMathJax) (u1), the choice-specific expected discounted values (StartMathJax)U_0(StopMathJax) (capU0) and (StartMathJax)U_1(StopMathJax) (capU1), and their contrast (StartMathJax)\Delta U(StopMathJax) (deltaU).

With (StartMathJax)\Delta U(StopMathJax) computed, and (StartMathJax)\Pi(StopMathJax) specified, we proceed to simulate a (StartMathJax)T\times N(StopMathJax) matrix of choices choices and a (StartMathJax)T\times N(StopMathJax) matrix of states iX (recall from Section 3 that iX contains indices that point to elements of (StartMathJax){\cal X}(StopMathJax) rather than those values themselves).

52Nested Fixed Point Maximum Likelihood Estimation

First, suppose that (StartMathJax)\Pi(StopMathJax) is known. We use fmincon from Matlab's Optimization Toolbox to maximize the partial likelihood for the choices (the code can easily be adapted to use other optimizers and packages, because these have a very similar syntax; see below). Because fmincon is a minimizer, we use minus the log likelihood as its objective. The function negLogLik computes this objective, but has input arguments other than the vector of model parameters to be estimated. Because the syntax of fmincon does not allow this, we define a function handle objectiveFunction to an anonymous function that equals negLogLik but does not have this extra inputs.

Before we can put fmincon to work on this objective function, we first have to set some of its other input arguments. We specify a (StartMathJax)3\times 1(StopMathJax) vector startvalues with starting values for the parameters to be estimated, (StartMathJax)(\beta_0,\beta_1,\delta_1)'(StopMathJax).

We also set a lower bound of 0 on the third parameter, (StartMathJax)\delta_1(StopMathJax), and (nonbinding) lower bounds of (StartMathJax)-\infty(StopMathJax) on the other two parameters (lowerBounds). There is no need to specify upper bounds.²

Finally, we pass some options, including tolerances that specify the criterion for the outer loop convergence, to fmincon through the structure OptimizerOptions (recall that we have already set the inner loop tolerance in tolFixedPoint). We use the function optimset from the Optimization Toolbox to assign values to specific fields (options) in OptimizerOptions and then call fmincon to run the NFXP maximum likelihood procedure (to use Knitro instead, simply replace fmincon by knitromatlab, knitrolink, or ktrlink, depending on the packages installed³).

This gives maximum partial likelihood estimates of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax). To calculate standard errors, we call negLogLik once more to estimate the corresponding Fisher information matrix and store this in informationMatrix. Its inverse is an estimate of the maximum likelihood estimator's asymptotic variance-covariance matrix.

The resulting parameter estimates and standard errors are displayed (third and fourth columns), together with the parameters' true (first column) and starting values (second column).

53Extension to an Unknown Markov Transition Matrix for the State

Finally, consider the more realistic case that (StartMathJax)\Pi(StopMathJax) is not known. In this case, Rust (1994) suggests a two-stage procedure. In the first stage, we estimate (StartMathJax)\Pi(StopMathJax) using estimatePi and store the results in a (StartMathJax)K\times K(StopMathJax) matrix piHat.

In the second stage, maximum partial likelihood estimates of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax) can be computed using the NFXP procedure, with piHat replacing capPi. This, with the question how the first-stage sampling error affects the precision of the second-stage estimator of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax), is left for the exercises.

61Theoretical Exercises

Theoretical Exercise 1

Prove that (StartMathJax)\Psi(StopMathJax) is a contraction.

Theoretical Exercise 2

Is (StartMathJax)\rho(StopMathJax) identified in the example model? What about (StartMathJax)\delta_0(StopMathJax)? (To the extent that these are identified, you may want to extend your numerical experiments below to include their estimation.) See Abbring (2010) for results and references.

Theoretical Exercise 3

Show the equivalence of policy iteration and the Newton-Kantorovich method for solving (StartMathJax)U-\Psi(U)=0(StopMathJax) for (StartMathJax)U(StopMathJax) (see Appendices 13 and 15). Read about the properties of policy iteration and discuss what these imply about the convergence properties of the Newton-Kantorovich method. Does policy iteration converge in a finite number of steps? Why (not)?

62Computational Exercises

The computational exercises typically ask you to modify the code so that it can handle alternative econometric procedures and models. This often requires that you adapt the analytic computation of the derivatives of the log likelihood. It may be convenient to first (or only) work with numerical derivatives by setting the option "GradObj" in fmincon or its alternative to "off"; this allows you to find the estimates without using the analytical gradient in the optimization step. Once you have coded up the analytical derivatives (in order to use them in the optimization step and to obtain an estimate of the standard errors), you are advised to verify them against numerical derivatives, for example by setting the option "DerivativeCheck" in fmincon to "on" (however, make sure to set this option to "off" once you have verified the derivatives and start using the code for repeated, e.g. Monte Carlo, calculations).

Computational Exercise 1

Extend the script in dynamicDiscreteChoice.m with a single simulation to a Monte Carlo experiment in which estimates are computed for a sequence of simulated data sets (for now, as in dynamicDiscreteChoice.m, take (StartMathJax)\Pi(StopMathJax) to be known). Keep track of the estimates and report statistics such as their means and standard deviations. Compare the latter to asymptotic standard errors. Make sure to set the seed of the random number generator so that your numerical results can be replicated. Note that you can also use this Monte Carlo setup in later questions to study the finite-sample performance of the various procedures studied.

Computational Exercise 2

For the case of estimating the demand for differentiated products using NFXP maximum likelihood, as in Berry, Levinsohn, and Pakes (1995), Dubé, Fox, and Su (2012) claim that it is important to carefully balance the convergence tolerances used in the inner (contraction) loop and the outer (likelihood maximization) loop. In particular, they argue that the inner loop needs to be computed at a much higher precision than the tolerance used in the outer loop. Experiment with the values of tolFixedPoint on the one hand and the tolerances in the tolX and tolFun fields of OptimizerOptions on the other hand to investigate this issue in the context of this package's firm entry and exit problem.

Computational Exercise 3

Implement the MPEC approach to maximum likelihood estimation of our structural model, as advocated by Judd and Su (2012).

• First, modify negLogLik so that it it takes (StartMathJax)(U_0,U_1)(StopMathJax) or (StartMathJax)\Delta U(StopMathJax) as an input argument (instead of solving the model for them) and computes the log (partial) likelihood directly from the choice probabilities implied by (StartMathJax)\Delta U(StopMathJax).
• Then, extend the script in dynamicDiscreteChoice.m so that it alternatively maximizes the log likelihood with respect to both the parameters of interest and the values of (StartMathJax)\Delta U(StopMathJax), subject to the constraint on (StartMathJax)\Delta U(StopMathJax) implied by (StartMathJax)U=\Psi(U)(StopMathJax) (which can be specified using the function bellman).

Would you expect the NFXP and MPEC approaches to give the same estimates of the parameters of interest (up to numerical precision)? How is the relative numerical performance of both procedures? Which one is faster? Compare your results to those in Iskhakov et al. (2016) and explain.

Computational Exercise 4

Compute the two-stage maximum partial likelihood estimates of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax) for the case that (StartMathJax)\Pi(StopMathJax) is not known. Are the estimators of the standard errors that assume (StartMathJax)\Pi(StopMathJax) known consistent if (StartMathJax)\Pi(StopMathJax) is estimated instead? Read Rust (1994) and explain why the two-stage estimator of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax) is not efficient. The full information maximum likelihood (FIML) estimator estimates all parameters, including the ones in (StartMathJax)\Pi(StopMathJax), simultaneously by maximizing the full likelihood for the observed state transitions and choices. Write an alternative for negLogLik that codes up this likelihood function and adapt the estimation script so that it obtains the FIML estimates and their estimated standard errors. Do not code up the gradient and do not use it in the likelihood maximization. Perform a Monte Carlo study of the two-stage and FIML estimators of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax). For both estimators, report the means and standard deviations of the estimates of (StartMathJax)(\beta_0,\beta_1,\delta_1)(StopMathJax) and the average estimates of the standard errors across Monte Carlo simulations. Compare and discuss. Finally, discuss the three-stage approach suggested by Rust (1994).

Computational Exercise 5

Implement a simulation-based two-step estimator along the lines of Hotz et al. (1994), Rust (1994), and Bajari, Benkard, and Levin (2007).

• Add a new function for nonparametrically estimating the choice probabilities (StartMathJax)p(a|X_t,A_{t-1})\equiv\Pr(A_t=a|X_t,A_{t-1})(StopMathJax) over the support of (StartMathJax)(X_t,A_{t-1})(StopMathJax). With enough data and a finite state space, you can simply use the conditional sample frequency.
• Add (to this same function or in a new function) a procedure for estimating (StartMathJax)\Delta U(StopMathJax) by inverting the estimated choice probabilities (as proposed by Hotz and Miller (1993) and Hotz et al. (1994)).
• Add a new function for simulating a given number of state (including (StartMathJax)\varepsilon_t(StopMathJax)) and choice paths of some given length from each state observed in the data and for each first choice that can be made at the beginning of the path. Make use of the linearity of the flow utility in the parameters to reduce the dimensionality of the objects returned by this function and needed for constructing the objective in the next step.
• Adapt negLogLik into an objective function that takes nonparametric estimates of (StartMathJax)p(StopMathJax) or (StartMathJax)\Delta U(StopMathJax) and (a relevant summary of) simulated states and choices as inputs. As possible objectives to be minimized, implement both a weighted distance between the nonparametric estimates of (StartMathJax)\Delta U(StopMathJax) and simulated values of (StartMathJax)\Delta U(StopMathJax), and minus a log pseudo-likelihood based on the choice probabilities implied by the simulated values of (StartMathJax)\Delta U(StopMathJax). When doing so exploit linearity of the flow utilities in the parameters.
• Finally, extend the script in dynamicDiscreteChoice.m so that it successively calls these functions to implement a two step procedure for estimating the model, as in Hotz and Miller (1993) and Hotz et al. (1994). Analyze the numerical and statistical performance of this procedure with both of the possible objective functions implemented and compare. Discuss the theoretical relation between both approaches; see Pesendorfer and Schmidt-Dengler (2008).

See the course slides for a brief and general description of this procedure and Rust (1994) for a detailed algorithm, with discussion. Note that the algorithm described in Rust (1994) is similar to that of Bajari, Benkard, and Levin (2007) for games that we will discuss later in the course. Both build on the ideas of Hotz et al. (1994), who use the special logit structure to simplify the simulation (see the discussion in Bajari, Benkard, and Levin (2007)).

Computational Exercise 6

Extend the code for simulation and NFXP estimation so that it can handle a finite number of unobserved types, as in the work of Eckstein and Wolpin (1990), Keane and Wolpin (1997), and Eckstein and Wolpin (1999). Specifically, suppose that there are two types in the population, one with entry cost (StartMathJax)\delta_1^1(StopMathJax) and the other with entry cost (StartMathJax)\delta_1^2(StopMathJax). The firms are distributed across both unobserved entry cost types independently from all other variables in the model. We would like to estimate (StartMathJax)\beta_0,\beta_1,\delta_1^1,\delta_1^2(StopMathJax), and the share of agents with entry cost (StartMathJax)\delta_1^1(StopMathJax). To this end:

• Extend simulateData so that it randomly draws an entry cost from a distribution with two points of support, (StartMathJax)\delta_1^1(StopMathJax) and (StartMathJax)\delta_1^2(StopMathJax), and simulates choice data corresponding to the entry cost drawn, for each firm, independently across firms.
• Extend negLogLik so that it computes each firm's likelihood contribution as the mixture (expectation) of the likelihood contributions for each entry cost type (which are given by the current specification of the likelihood for homogeneous firms) over the distribution of these two types in the population of firms.
• Finally, extend the script in dynamicDiscreteChoice.m so that it successively calls these functions. Check whether you can recover the entry cost distribution, and the other parameters, well.

Computational Exercise 7

Implement a version of the two-step estimator of Hotz et al. (1994) that similarly allows for a finite number of unobserved types. To this end, combine this estimator with the EM algorithm as in Arcidiacono and Miller (2011) (see also Arcidiacono and Ellickson (2011)).

Computational Exercise 8

Perform a Monte Carlo study of all the estimators that you have implemented. Evaluate and compare their numerical performance as measured in terms of convergence success and computation time; and their finite sample statistical performance in terms of, among other things, mean squared errors. Experiment with different choices of the numerical design parameters, such as the inner and outer loop tolerances used in the NFXP procedure.

Computational Exercise 9

Extend fixedPoint so that it has the option to compute (StartMathJax)U(StopMathJax) using the Newton-Kantorovich method instead of successive approximations, or a hybrid method that combines both (as in Rust (1987)). Note that, with a finite state space as in our example, the Newton-Kantovorich method reduces to the Newton-Raphson method and requires the Jacobian (StartMathJax)\partial\tilde\Psi_\theta(\bar U)/\partial\bar U'(StopMathJax) computed in an intermediate step of negLogLik. Explore the numerical performance of the various methods. What is your preferred method?

• Abbring, Jaap H. (2010). Identification of dynamic discrete choice models. Annual Review of Economics 2, 367-394.
• Aguirregabiria, Victor and Pedro Mira (2002). Swapping the nested fixed point algorithm: A class of estimators for discrete Markov decision models. Econometrica 70(4), 1519-1543.
• Arcidiacono, Peter and Paul B. Ellickson (2011). Practical methods for estimation of dynamic discrete choice models. Annual Review of Economics 3, 363-394.
• Arcidiacono, Peter and Robert A. Miller (2011). Conditional choice probability estimation of dynamic discrete choice models with unobserved heterogeneity. Econometrica 79, 1823-1867.
• Bajari, Patrick, Lanier Benkard, and Jonathan Levin (2007). Estimating dynamic models of imperfect competition. Econometrica 75(5), 1331-1370.
• Berry, Steven, James Levinsohn, and Ariel Pakes (1995). Automobile prices in market equilibrium. Econometrica 63, 841-890.
• Dixit, Avinash K. (1989). Entry and exit decisions under uncertainty. Journal of Political Economy 97(3), 620-638.
• Dubé, Jean-Pierre, Jeremy T. Fox, and Che-Lin Su (2012). Improving the numerical performance of static and dynamic aggregate discrete choice random coefficients demand estimation. Econometrica 80, 2231-2267.
• Eckstein, Zvi and Kenneth I. Wolpin (1990). Estimating a market equilibrium search model from panel data on individuals. Econometrica 58(4), 783-808.
• Eckstein, Zvi and Kenneth I. Wolpin (1999). Why youths drop out of high school: The impact of preferences, opportunities, and abilities. Econometrica 67(6), 1295-1339.
• Hotz, V. Joseph and Robert A. Miller (1993). Conditional choice probabilities and the estimation of dynamic models. Review of Economic Studies 60, 497-529.
• Hotz, V. Joseph, Robert A. Miller, Seth Sanders and Jeffrey Smith (1994). A simulation estimator for dynamic models of discrete choice. Review of Economic Studies 61(2), 265-289.
• Iskhakov, Fedor, Jinhyuk Lee, John Rust, Bertel Schjerning and Kyoungwon Seo (2016). Comment on "Constrained optimization approaches to estimation of structural models". Econometrica 84(1), 365-370.
• Judd, Kenneth L. 1998. Numerical Methods in Economics. MIT Press. Cambridge, MA.
• Judd, Kenneth L. and Che-Lin Su (2012). Constrained optimization approaches to estimation of structural models. Econometrica 80(5), 2213-2230.
• Keane, Michael P. and Kenneth I. Wolpin (1997). The career decisions of young men. Journal of Political Economy 105(3), 473-522.
• Ljungqvist, Lars and Thomas J. Sargent (2000). Recursive Macroeconomic Theory, Second edition. MIT Press. Cambridge, MA.
• McFadden, Daniel (1981). Econometric models of probabilistic choice. In C. Manski and D. McFadden (Eds.). Discrete Data with Econometric Applications. MIT Press. Cambridge, MA.
• Pesendorfer, Martin and Philipp Schmidt-Dengler (2008). Asymptotic least squares estimators for dynamic games. Review of Economic Studies 75(3), 901-928.
• Puterman, Martin L. and Shelby L. Brumelle (1979). On the convergence of policy iteration in stationary dynamic programming. Mathematics of Operations Research 4.1, 60-69.
• Rust, John (1987). Optimal replacement of GMC bus engines: An empirical model of Harold Zurcher. Econometrica 55, 999-1033.
• Rust, John (1994). Structural estimation of Markov decision processes. In R. Engle and D. McFadden (Eds.). Handbook of Econometrics 4, 3081-3143. North-Holland. Amsterdam.
• Stokey, Nancy L. and Robert E. Lucas (with Edward C. Prescott) (1989). Recursive Methods in Economic Dynamics. Harvard University Press. Cambridge, MA.

11Definitions

A metric space is a set (StartMathJax)\cal{U}(StopMathJax) with a metric (StartMathJax)m:\cal{U}\times\cal{U}\rightarrow\mathbb{R}(StopMathJax) such that, for any (StartMathJax)U,U',U''\in\cal{U}(StopMathJax),

1. (StartMathJax)m(U,U')=0(StopMathJax) if and only if (StartMathJax)U=U'(StopMathJax),
2. (StartMathJax)m(U,U')=m(U',U)(StopMathJax), and
3. (StartMathJax)m(U,U'')\leq m(U,U')+m(U',U'')(StopMathJax) (triangle inequality).

A Cauchy sequence in (StartMathJax)({\cal U},m)(StopMathJax) is a sequence (StartMathJax)\{U_n\}(StopMathJax) in (StartMathJax){\cal U}(StopMathJax) such that, for each (StartMathJax)\epsilon \gt 0(StopMathJax), there exists some (StartMathJax)n(\epsilon)\in\mathbb{N}(StopMathJax) such that (StartMathJax)m(U_n,U_{n'}) \lt \epsilon(StopMathJax) for all (StartMathJax)n,n'\geq n(\epsilon)(StopMathJax).

A metric space (StartMathJax)({\cal U},m)(StopMathJax) is complete if every Cauchy sequence in (StartMathJax)({\cal U},m)(StopMathJax) converges to a point in (StartMathJax){\cal U}(StopMathJax).

A map (StartMathJax)\Psi:{\cal U}\rightarrow{\cal U}(StopMathJax) is a contraction with modulus (StartMathJax)\rho\in(0,1)(StopMathJax) if (StartMathJax)m\left[\Psi(U),\Psi(U')\right]\leq \rho m(U,U')(StopMathJax) for all (StartMathJax)U,U'\in{\cal U}(StopMathJax).

12Contraction Mapping (Banach Fixed Point) Theorem

If (StartMathJax)({\cal U},m)(StopMathJax) is a complete metric space and (StartMathJax)\Psi:{\cal U}\rightarrow{\cal U}(StopMathJax) is a contraction, then there exists a unique (StartMathJax)U\in{\cal U}(StopMathJax) such that (StartMathJax)U=\Psi(U)(StopMathJax). Furthermore, for any (StartMathJax)U_0\in{\cal U}(StopMathJax), the sequence (StartMathJax)\{U_n\}(StopMathJax) with (StartMathJax)U_{n+1}=\Psi(U_{n})(StopMathJax), (StartMathJax)n\in\mathbb{N}(StopMathJax), converges to (StartMathJax)U(StopMathJax).

Sketch of Proof

Suppose (StartMathJax)U,U'\in{\cal U}(StopMathJax) are such that (StartMathJax)U=\Psi(U)(StopMathJax) and (StartMathJax)U'=\Psi(U')(StopMathJax). Then, (StartMathJax)0\leq m(U,U')=m\left[\Psi(U),\Psi(U')\right]\leq\rho m(U,U')(StopMathJax) for some (StartMathJax)\rho\in(0,1)(StopMathJax). Consequently, (StartMathJax)m(U,U')=0(StopMathJax) and (StartMathJax)U=U'(StopMathJax).

Because (StartMathJax)\Psi(StopMathJax) is a contraction, any (StartMathJax)\{U_n\}(StopMathJax) it generates is a Cauchy sequence. Because (StartMathJax)({\cal U},m)(StopMathJax) is complete, it converges to some (StartMathJax)U\in{\cal U}(StopMathJax). Because contractions are continuous, (StartMathJax)U=\Psi(U)(StopMathJax).

13Computing the Fixed Point of a Contraction

The method of successive approximations directly applies the Contraction Mapping Theorem and approximates the fixed point (StartMathJax)U(StopMathJax) of a contraction (StartMathJax)\Psi(StopMathJax) with the sequence (StartMathJax)\{U_n\}(StopMathJax) generated from (StartMathJax)U_{n+1}=\Psi(U_n)(StopMathJax) on a finitely discretized state space. This method has global but linear convergence.

The Newton-Kantorovich method searches for a zero of (StartMathJax)I-\Psi(StopMathJax), where (StartMathJax)I:U\in{\cal U}\mapsto U(StopMathJax) is the identity mapping. It approximates (StartMathJax)U(StopMathJax) with the sequence generated from

with (StartMathJax)I-\Psi'_{U_n}(StopMathJax) the Fréchet derivative of (StartMathJax)I-\Psi(StopMathJax) at (StartMathJax)U_n(StopMathJax) (with a discrete state space, this is simply the linear map defined by the appropriate finite dimensional Jacobian, as in negLogLik).

14Blackwell's Sufficient Conditions for a Contraction

Let (StartMathJax){\cal U}(StopMathJax) be the space of functions (StartMathJax)U:{\cal X}\rightarrow\mathbb{R}(StopMathJax) such that (StartMathJax)\sup|U| \lt \infty(StopMathJax) ((StartMathJax)U(StopMathJax) bounded), with metric (StartMathJax)m(U,U')=\sup|U-U'|(StopMathJax). Suppose that, for some (StartMathJax)\rho\in(0,1)(StopMathJax), (StartMathJax)\Psi:{\cal U}\rightarrow{\cal U}(StopMathJax) satisfies

1. (monotonicity) (StartMathJax)\Psi(U)\leq \Psi(U')(StopMathJax) and
2. (discounting) (StartMathJax)\Psi(U+\gamma)\leq \Psi(U)+\rho\gamma(StopMathJax)

for all (StartMathJax)U,U'\in{\cal U}(StopMathJax) such that (StartMathJax)U\leq U'(StopMathJax) and all (StartMathJax)\gamma\in(0,\infty)(StopMathJax). Then, (StartMathJax)\Psi(StopMathJax) is a contraction with modulus (StartMathJax)\rho(StopMathJax).

Sketch of Proof

Let (StartMathJax)U,U'\in{\cal U}(StopMathJax) be such that (StartMathJax)U\leq U'(StopMathJax). Note that (StartMathJax)U\leq U'+m(U,U')(StopMathJax), so that (StartMathJax)\Psi(U)\leq \Psi(U')+\rho m(U,U')(StopMathJax) by monotonicity and discounting. Similarly, (StartMathJax)\Psi(U')\leq \Psi(U)+\rho m(U,U')(StopMathJax). Thus, (StartMathJax)m\left[\Psi(U),\Psi(U')\right]\leq\rho m(U,U')(StopMathJax).

15Application to Dynamic Programming

The choice-specific value function (StartMathJax)U(StopMathJax) is a fixed point of (StartMathJax)\Psi(StopMathJax), which satisfies Blackwell's sufficient conditions with (StartMathJax)\rho(StopMathJax) equal to the discount factor.

• Under regularity conditions such that (StartMathJax)U(StopMathJax) is bounded, this ensures that it is the unique fixed point of a contraction on the complete (Banach) space (StartMathJax)({\cal U},m)(StopMathJax) of bounded functions.
• If (StartMathJax){\cal U}'\subset{\cal U}(StopMathJax) is closed and (StartMathJax)\Psi(U')\in{\cal U}'(StopMathJax) for all (StartMathJax)U'\in{\cal U}'(StopMathJax), then (StartMathJax)\Psi(StopMathJax) is a contraction on the complete subspace (StartMathJax)({\cal U}',m)(StopMathJax), so that (StartMathJax)U\in{\cal U}'(StopMathJax) (examples: monotonicity, continuity, etcetera).
• (StartMathJax)U(StopMathJax) can be computed by successive approximations, the Newton-Kantorovich method, or a hybrid algorithm as in Rust (1987).

In this context, the method of successive approximations is alternatively referred to as value iteration. Moreover, the Newton-Kantorovich method is closely related to an alternative method called policy iteration. Each of this method's iterations takes a value function as input, determines the corresponding optimal policy, and then gives the values of applying that policy forever as output (see e.g. Ljungqvist and Sargent (2000), Section 3.1.1, which refers to this method as Howard's policy improvement algorithm, or Rust (1994), Section 2.5).

Puterman and Brumelle (1979) show that policy iteration is generally equivalent to applying the Newton-Kantorovich method to finding the fixed point of the Bellman equation (in the sense that both produce the same sequence of value functions when starting from the same value). Ljungqvist and Sargent (2000), Section 4.4, present and develop this result for the special case of a finite state space.

Finally, for our example's special case of a finite (StartMathJax){\cal X}(StopMathJax), Aguirregabiria and Mira (2002)'s Proposition 1 establishes that policy iteration is equivalent to applying the Newton-Kantorovich method to finding a fixed point of the Bellman-like operator (StartMathJax)\Psi(StopMathJax) (which Aguirregabiria and Mira (2002) refer to as a smoothed Bellman operator).