Dynamic Discrete Choice Models: Methods, Matlab Code, and Exercises

Jaap H. Abbring§Tobias J. Klein

November 2023

This document supports the first computing sessions in a graduate course on dynamic discrete choice models. It is centered around basic Matlab code for solving, simulating, and empirically analyzing a simple model of firm entry and exit under uncertainty. This code is available from a public Github repository and can be distributed to students as starter code, for example using the Github Classroom. Exercises ask students to adapt and extend this starter code to apply different and more advanced computational and econometric methods to a wider range of models.

We thank Jeffrey Campbell for help with Komments++ and Nan Yang and our students in CentER Tilburg's Empirical Industrial Organization 2 (230323) for comments on earlier versions of this code. First draft: February 8, 2011. We welcome the use of this software under an MIT license. DOI: 10.5281/zenodo.4287733. © 2020 Jaap H. Abbring and Tobias J. Klein.
§Department of Econometrics & OR, Tilburg University, P.O. Box 90153, 5000 LE Tilburg, The Netherlands. E-mail: jaap@abbring.org. Web: jaap.abbring.org.
Department of Econometrics & OR, Tilburg University, P.O. Box 90153, 5000 LE Tilburg, The Netherlands. E-mail: t.j.klein@uvt.nl. Web: www.tobiasklein.ws.

Viewing and Using this File

This file (dynamicDiscreteChoice.m.html) documents the code in the GitHub repository jabbring/dynamic-discrete-choice (a Zip archive packages both together). It was generated from the Matlab script dynamicDiscreteChoice.m using the Komments++ package, which was created and generously provided to us by Jeffrey R. Campbell. It documents how you can run the script dynamicDiscreteChoice.m with Matlab to specify, simulate, and estimate an empirical discrete-time version of the model of firm entry and exit under uncertainty by Dixit (1989). These computations will use, and therefore exemplify the use of, various tailor-made Matlab functions that are documented in this file. Thus, this file also is a guide to these functions and the way they can be adapted and used in your own exercises.

In Safari and Firefox, you can switch between the default view of this document, which displays the working code with all its documentation, and an alternative view that shows only the code by pressing “c”.

Software Requirements

Running the code documented in this file requires a recent version of Matlab with its Optimization Toolbox. The code can easily be adapted to use the Knitro solver instead of Matlab's Optimization Toolbox (see Section 5).1

Implementations in Other Languages

A Julia implementation of this code by Rafael Greminger is available from the GitHub repository rgreminger/DDCModelsExample.jl.

Road Map

The remainder of this document proceeds as follows. The next section covers two functions that define the decision problem, flowpayoffs and bellman. The versions of these functions handed out to the students, and documented here, define a very basic entry and exit problem with sunk costs and ongoing uncertainty. By changing these functions, different and more extensive decision problems can be studied with the procedures documented in later sections. Section 2 discusses fixedPoint, which solves for the optimal choice-specific expected discounted values (and therewith for the optimal decision rule) by finding the unique fixed point of the contraction mapping defined by bellman. Section 3 documents how simulateData can be used to simulate data for the decision problem set up by flowpayoffs and bellman. Section 4 documents functions for estimating the transition matrix (StartMathJax)\Pi(StopMathJax) and for computing the log partial likelihood for conditional choice data. Section 5 documents the script that uses these functions to set up the model, simulate data, and estimate the model with Rust (1987)'s nested fixed point (NFXP) maximum likelihood method. Section 6 contains a range of student exercises. Appendix 1 provides mathematical background by discussing some theory of contractions. Appendix 2 documents procedures that are called by the main routines, but that are of limited substantial interest, such as randomDiscrete.

1A Firm Entry and Exit Problem

Consider a simple, discrete-time version of the model of firm entry and exit under uncertainty in Dixit (1989). At each time (StartMathJax)t\in\mathbb{N}(StopMathJax), a firm can choose to either serve a market (choice 1) or not (choice 0). The payoffs from either choice depend on the firm's choice (StartMathJax)A_{t-1}(StopMathJax) in the previous period, because the firm may have to incur a sunk cost to enter or exit the market (for definiteness, suppose that the firm is initially inactive, (StartMathJax)A_0=0(StopMathJax)). Payoffs also depend on externally determined (choice-independent) observed scalar state variables (StartMathJax)X_t(StopMathJax), as well as unobserved state variables (StartMathJax)\varepsilon_{t}(0)(StopMathJax) and (StartMathJax)\varepsilon_{t}(1)(StopMathJax) that are revealed to the firm before it makes its choice in period (StartMathJax)t(StopMathJax). Specifically, in period (StartMathJax)t(StopMathJax), its flow profits from choice 0 are

(0)

(StartMathJax)u_0(X_t,A_{t-1})+\varepsilon_{t}(0)= -A_{t-1}\delta_0+\varepsilon_{t}(0),(StopMathJax)

and its flow profits from choice 1 are

(0)

(StartMathJax)u_1(X_t,A_{t-1})+\varepsilon_{t}(1)=\beta_0+\beta_1 X_t-(1-A_{t-1})\delta_1+\varepsilon_{t}(1).(StopMathJax)

Note that (StartMathJax)\delta_0(StopMathJax) and (StartMathJax)\delta_1(StopMathJax) are exit and entry costs that the firm only pays if it changes state. Gross of these costs, an active firm makes profits (StartMathJax)\beta_0+\beta_1 X_t+\varepsilon_{t}(1)(StopMathJax) and an inactive firm has profits (StartMathJax)\varepsilon_{t}(0)(StopMathJax).

The state variable (StartMathJax)X_t(StopMathJax) has finite support (StartMathJax){\cal X}\equiv\{x^1,\ldots,x^K\}(StopMathJax). From its random initial value (StartMathJax)X_0(StopMathJax), it follows a first-order Markov chain with (StartMathJax)K\times K(StopMathJax) transition probability matrix (StartMathJax)\Pi(StopMathJax), with typical element (StartMathJax)\Pi_{ij}=\Pr(X_{t+1}=x^j|X_t=x^i)(StopMathJax), independently of the firm's choices. The profit shocks (StartMathJax)\varepsilon_{t}(a)(StopMathJax) are independent of (StartMathJax)\{X_t\}(StopMathJax) and across time (StartMathJax)t(StopMathJax) and choices (StartMathJax)a(StopMathJax), with type I extreme value distributions centered around 0. Like (StartMathJax)X_t(StopMathJax), they may affect but are not affected by the firm's choices. Thus, (StartMathJax)X_t(StopMathJax) and (StartMathJax)\varepsilon_t\equiv[\varepsilon_t(0),\varepsilon_t(1)](StopMathJax) are “exogenous” state variables of which the evolution is externally specified. Consequently, the firm controls the evolution of the state (StartMathJax)(X_t,A_{t-1},\varepsilon_t)(StopMathJax) only through its choice of (StartMathJax)A_{t-1}(StopMathJax).

The firm has rational expectations about future states and choices. It chooses the action (StartMathJax)A_t(StopMathJax) that maximizes its expected flow of profits, discounted at a factor (StartMathJax)\rho \lt 1(StopMathJax).

Two functions, flowpayoffs and bellman, together code up this simple model. If you wish to experiment with other functional forms for the flow profits, you should edit flowpayoffs.m. The model's dynamic specification can be changed in bellman.m.

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).

flowpayoffs.m7function [u0,u1] = flowpayoffs(supportX,beta,delta)

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.

flowpayoffs.m24nSuppX = size(supportX,1);

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

(0)

(StartMathJax) \left[\begin{array}{ccc} u_0(x^1,0)&~~~&u_0(x^1,1)\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ u_0(x^{K},0)&~~~&u_0(x^{K},1) \end{array}\right] = \left[\begin{array}{ccc} 0&~~~&-\delta_0\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ 0&~~~&-\delta_0 \end{array}\right] (StopMathJax)

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

(0)

(StartMathJax) \left[\begin{array}{ccc} u_1(x^1,0)&~~~&u_1(x^1,1)\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ u_1(x^K,0)&~~~&u_1(x^{K},1) \end{array}\right] = \left[\begin{array}{ccc} \beta_0+\beta_1 x^1-\delta_1&~~~&\beta_0+\beta_1 x^1\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ \cdot&~~~&\cdot\\ \beta_0+\beta_1 x^{K}-\delta_1&~~~&\beta_0+\beta_1 x^{K} \end{array}\right]. (StopMathJax)

flowpayoffs.m55u0 = [zeros(nSuppX,1) -delta(1)*ones(nSuppX,1)];
flowpayoffs.m56u1 = [ones(nSuppX,1) supportX]*beta*[1 1]-delta(2)*ones(nSuppX,1)*[1 0];

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

(1)

(StartMathJax)\label{eq:bellman} U_a(X_t,A_{t-1})=u_a(X_t,A_{t-1})+\rho\mathbb{E}\left[R_a(X_{t+1})|X_t\right], (StopMathJax)

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

(2)

(StartMathJax)\label{eq:surplus} R_a(x)=\log\left\{\exp\left[U_0(x,a)\right]+\exp\left[U_1(x,a)\right]\right\}. (StopMathJax)

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

bellman.m15function [capU0,capU1] = bellman(capU0,capU1,u0,u1,capPi,rho)

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.

bellman.m34r0 = log(exp(capU0(:,1))+exp(capU1(:,1)));
bellman.m35r1 = log(exp(capU0(:,2))+exp(capU1(:,2)));

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

bellman.m39capU0 = u0 + rho*capPi*r0*[1 1];
bellman.m40capU1 = u1 + rho*capPi*r1*[1 1];

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).

2Solving the Decision Problem

The function fixedPoint computes the fixed point (StartMathJax)U(StopMathJax) of the Bellman-like operator (StartMathJax)\Psi(StopMathJax), using the method of successive approximations. It is easy to show that (StartMathJax)\Psi(StopMathJax) is a contraction, so that the method of successive approximations converges linearly and globally to a point within a positive maximum absolute distance tolFixedPoint (see Appendix 13).

fixedPoint.m7function [capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,bellman,capU0,capU1)

It requires the following input arguments:

{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);
{rho} a scalar with the value of the discount factor (StartMathJax)\rho(StopMathJax);
{tolFixedPoint} a scalar tolerance level that is used to determine convergence of the successive approximations;
{bellman} the handle to the function [capU0,capU1] = bellman(capU0,capU1,u0,u1,capPi,rho) that iterates once on (StartMathJax)\Psi(StopMathJax);
{capU0} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is a starting value for (StartMathJax)U_0(x^i,j-1)(StopMathJax) (optional; set to [] to select default starting value);
{capU1} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is a starting value for (StartMathJax)U_1(x^i,j-1)(StopMathJax) (optional; set to [] to select default starting value).

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); and
{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).

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

fixedPoint.m28nSuppX = size(capPi,1);

The starting values for (StartMathJax)U_0(StopMathJax) and (StartMathJax)U_1(StopMathJax) are set to (StartMathJax)0(StopMathJax) if the input arguments capU0 and capU1 are empty.

fixedPoint.m32if isempty(capU0)
fixedPoint.m33 capU0 = zeros(nSuppX,2);
fixedPoint.m34end
fixedPoint.m35if isempty(capU1)
fixedPoint.m36 capU1 = zeros(nSuppX,2);
fixedPoint.m37end

The (StartMathJax)K\times 2(StopMathJax) matrices inU0 and inU1 store the values of (StartMathJax)U(StopMathJax) that are fed into the operator (StartMathJax)\Psi(StopMathJax), for comparison with the value of (StartMathJax)\Psi(U)(StopMathJax) that is subsequently stored in capU0 and capU1. They are initialized to deviate from capU0 and capU1 by more than tolFixedPoint, so that the while statement allows at least one iteration of (StartMathJax)\Psi(StopMathJax), and stops as soon as (StartMathJax)\max\{\max|\Psi_0(U)-U_0|,\max|\Psi_1(U)-U_1|\}(StopMathJax) no longer exceeds the tolerance level in tolFixedPoint.

fixedPoint.m41inU0 = capU0+2*tolFixedPoint;
fixedPoint.m42inU1 = capU1+2*tolFixedPoint;
fixedPoint.m43while (max(max(abs(inU0-capU0)))>tolFixedPoint) || (max(max(abs(inU1-capU1)))>tolFixedPoint);
fixedPoint.m44 inU0 = capU0;
fixedPoint.m45 inU1 = capU1;
fixedPoint.m46 [capU0,capU1] = bellman(inU0,inU1,u0,u1,capPi,rho);
fixedPoint.m47end

You can replace fixedPoint by another function if you want to use alternative methods, such as Newton methods, for computing the fixed point (StartMathJax)U(StopMathJax), or if you want to work with finite horizon problems.

3Data Simulation

Suppose that we have computed (StartMathJax)\Delta U(x,a)\equiv U_1(x,a)-U_0(x,a)(StopMathJax) for all (StartMathJax)(x,a)\in{\cal X}\times\{0,1\}(StopMathJax), using e.g. flowpayoffs and fixedPoint, and that we have specified the Markov transition matrix (StartMathJax)\Pi(StopMathJax). Then, the function simulateData can be used to simulate (StartMathJax)N(StopMathJax) independent histories (StartMathJax)\{(X_1,A_1),\ldots,(X_T,A_T)\}(StopMathJax), taking (StartMathJax)A_0=0(StopMathJax) as given and drawing (StartMathJax)X_1(StopMathJax) from the stationary distribution of (StartMathJax)\{X_t\}(StopMathJax).

simulateData.m7function [choices,iX] = simulateData(deltaU,capPi,nPeriods,nFirms)

The function simulateData requires the following input arguments:

{deltaU} a (StartMathJax)K\times 2(StopMathJax) matrix of which the (StartMathJax)(i,j)(StopMathJax)th entry is (StartMathJax)\Delta U(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);
{nPeriods} the scalar number (StartMathJax)T(StopMathJax) of time periods to simulate data for; and
{nFirms} the scalar number (StartMathJax)N(StopMathJax) of firms to simulate data for.

It returns

{choices} a (StartMathJax)T\times N(StopMathJax) matrix with simulated choices, with each each column containing an independent simulation of (StartMathJax)(A_1,\ldots,A_T)(StopMathJax); and
{iX} a (StartMathJax)T\times N(StopMathJax) matrix with simulated states, with each each column containing the indices (in (StartMathJax){\cal X}(StopMathJax)) of an independent simulation of (StartMathJax)(X_1,\ldots,X_T)(StopMathJax). For example, if (StartMathJax)X_1=x^3(StopMathJax) for some firm, then 3, not (StartMathJax)x^3(StopMathJax), is stored in iX.

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

simulateData.m23nSuppX = size(capPi,1);

Next, it assumes that (StartMathJax)\{X_t\}(StopMathJax) is ergodic and initializes the simulation by drawing (StartMathJax)X_1(StopMathJax), (StartMathJax)N(StopMathJax) independent times, from the stationary distribution of (StartMathJax)\{X_t\}(StopMathJax). To this end, it first solves

(0)

(StartMathJax)\left[\begin{array}{cccc} 1-\Pi_{11} &-\Pi_{21} &\cdots &-\Pi_{K1}\\ -\Pi_{12} &1-\Pi_{22}&\cdots &-\Pi_{K2}\\ \vdots & &\ddots &\vdots\\ -\Pi_{1(K-1)}&\cdots &~~1-\Pi_{(K-1)(K-1)}~~&-\Pi_{K(K-1)}\\ 1 &\cdots &1 &1 \end{array}\right]P^\infty= \left(\begin{array}{c} 0\\0\\ \vdots\\ 0\\1 \end{array}\right) (StopMathJax)

for the (StartMathJax)K\times 1(StopMathJax) vector (StartMathJax)P^\infty(StopMathJax) with the stationary probabilities (StartMathJax)P^\infty_k=\lim_{t\rightarrow\infty}\Pr(X_t=x^k)(StopMathJax) and stores the result in pInf.

simulateData.m39oneMinPi = eye(nSuppX)-capPi';
simulateData.m40pInf = [oneMinPi(1:nSuppX-1,:);ones(1,nSuppX)]\[zeros(nSuppX-1,1);1];

Then, it uses the auxiliary function randomDiscrete (see Appendix 2) and the values stored in pInf to simulate a (StartMathJax)1\times N(StopMathJax) vector of values of (StartMathJax)X_1(StopMathJax) from the stationary distribution (StartMathJax)P^\infty(StopMathJax) and stores their indices in iX.

simulateData.m44iX = randomDiscrete(pInf*ones(1,nFirms));

Using these (StartMathJax)N(StopMathJax) simulated values of (StartMathJax)X_1(StopMathJax), and (StartMathJax)N(StopMathJax) simulated values of (StartMathJax)-\Delta\varepsilon_1\equiv\varepsilon_1(0)-\varepsilon_1(1)(StopMathJax) that are stored in deltaEpsilon, it simulates (StartMathJax)N(StopMathJax) values of the first choice by using that (StartMathJax)A_1=1(StopMathJax) if (StartMathJax)\Delta U(X_1,0) \gt -\Delta\varepsilon_1(StopMathJax) and (StartMathJax)A_1=0(StopMathJax) otherwise. These are stored in the (StartMathJax)1\times N(StopMathJax) vector choices.

simulateData.m48deltaEpsilon = random('ev',zeros(1,nFirms),ones(1,nFirms))-random('ev',zeros(1,nFirms),ones(1,nFirms));
simulateData.m49choices = deltaU(iX,1)' > deltaEpsilon;

Finally, (StartMathJax)N(StopMathJax) values of (StartMathJax)X_t(StopMathJax) are simulated, using the transition matrix (StartMathJax)\Pi(StopMathJax) and randomDiscrete, and their indices added as a row to the bottom of the (StartMathJax)(t-1)\times N(StopMathJax) matrix iX; and (StartMathJax)N(StopMathJax) values of (StartMathJax)A_t(StopMathJax) are simulated, using that (StartMathJax)A_t=1(StopMathJax) if (StartMathJax)\Delta U(X_t,A_{t-1}) \gt -\Delta\varepsilon_t(StopMathJax) and (StartMathJax)A_t=0(StopMathJax) otherwise, and stored as a row at the bottom of the (StartMathJax)(t-1)\times N(StopMathJax) matrix choices; recursively for (StartMathJax)t=2,\ldots,T(StopMathJax).

simulateData.m53for t = 2:nPeriods
simulateData.m54 iX = [iX;randomDiscrete(capPi(iX(end,:),:)')];
simulateData.m55 deltaEpsilon = random('ev',zeros(1,nFirms),ones(1,nFirms))-random('ev',zeros(1,nFirms),ones(1,nFirms));
simulateData.m56 choices = [choices;(deltaU(iX(end,:)+nSuppX*choices(end,:)) > deltaEpsilon)];
simulateData.m57end

4Estimation

Suppose that we have a random sample (StartMathJax)\left\{\left[(x_{1n},a_{1n}),\ldots,(x_{Tn},a_{Tn})\right]; n=1,\ldots,N\right\}(StopMathJax) from the population of state and choice histories (StartMathJax)\{(X_1,A_1),\ldots,(X_T,A_T)\}(StopMathJax). Following Rust (1994), we provide functions for implementing a two-stage procedure in which, first, the state transition matrix (StartMathJax)\Pi(StopMathJax) is estimated directly from observed state transitions and, second, the remaining parameters are estimated by maximizing the partial likelihood for the observed choices.

41First Stage: State Transitions

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

estimatePi.m7function piHat = estimatePi(iX,nSuppX)

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.

estimatePi.m20nPeriods=size(iX,1);

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.

estimatePi.m24for i=1:nSuppX
estimatePi.m25 for j=1:nSuppX
estimatePi.m26 piHat(i,j)=sum(sum((iX(2:nPeriods,:)==j)&(iX(1:nPeriods-1,:)==i)))/sum(sum((iX(1:nPeriods-1,:)==i)));
estimatePi.m27 end
estimatePi.m28end

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).

negLogLik.m7function [nll,negScore,informationMatrix] =
negLogLik.m8 negLogLik(choices,iX,supportX,capPi,beta,delta,rho,flowpayoffs,bellman,fixedPoint,tolFixedPoint)

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.

negLogLik.m35nSuppX = size(supportX,1);

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.

negLogLik.m39[u0,u1] = flowpayoffs(supportX,beta,delta);
negLogLik.m40[capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,bellman,[],[]);
negLogLik.m41deltaU = capU1-capU0;
negLogLik.m42pExit = 1./(1+exp(deltaU));

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

(0)

(StartMathJax)p(a_{tn}|x_{tn},a_{(t-1)n})=a_{tn}+\frac{1-2 a_{tn} }{1+\exp\left[\Delta U(x_{tn},a_{(t-1)n})\right]},(StopMathJax)

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.

negLogLik.m49laggedChoices = [zeros(1,size(choices,2));choices(1:end-1,:)];
negLogLik.m50p = choices + (1-2*choices).*pExit(iX+nSuppX*laggedChoices);
negLogLik.m51nll = -sum(sum(log(p)));

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.

negLogLik.m57if nargout>=2

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

(0)

(StartMathJax)\frac{\partial\log\left[\prod_{t=1}^Tp(a_{tn}|x_{tn},a_{(t-1)n})\right]}{\partial \theta} = -\sum_{t=1}^T\left(1-2 a_{tn}\right)\left[1-p(a_{tn}|x_{tn},a_{(t-1)n})\right] \frac{\partial\Delta U(x_{tn},a_{(t-1)n})}{\partial\theta}. (StopMathJax)

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,

(0)

(StartMathJax) \bar U=\left[U_0(x^1,0),\ldots,U_0(x^K,0),U_0(x^1,1),\ldots,U_0(x^K,1),U_1(x^1,0),\ldots,U_1(x^K,0),U_1(x^1,1),\ldots,U_1(x^K,1)\right]', (StopMathJax)

and that satisfies

(0)

(StartMathJax)\bar U= \tilde\Psi_\theta(\bar U),(StopMathJax)

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

(0)

(StartMathJax)\left[I_{4K} - \frac{\partial \tilde \Psi_\theta(\bar U)}{\partial\bar U'}\right]\frac{\partial\bar U}{\partial \theta'} = \frac{\partial \tilde\Psi_\theta(\bar U)}{\partial\theta'}(StopMathJax)

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;

(0)

(StartMathJax)\frac{\partial \tilde\Psi_\theta(\bar U)}{\partial\bar U'}= \left(\begin{array}{cccc} ~D_{00}~ &~O_K~ &~D_{10}~ &~O_K~\\ ~D_{00}~ &~O_K~ &~D_{10}~ &~O_K~\\ ~O_K~ &~D_{01}~ &~O_K~ &~D_{11}~\\ ~O_K~ &~D_{01}~ &~O_K~ &~D_{11}~ \end{array}\right), (StopMathJax)

with

(0)

(StartMathJax)D_{a'a}=\rho\Pi\left[\begin{array}{cccc} p(a'|x^1,a) &0 &\cdots &0\\ 0 &p(a'|x^2,a)& &\vdots\\ \vdots & &\ddots &0\\ 0 &\cdots &0 &p(a'|x^K,a) \end{array}\right], (StopMathJax)

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

(0)

(StartMathJax)\frac{\partial \tilde \Psi_\theta(\bar U)}{\partial\theta'}= \left[ \begin{array}{cccc} ~0~ & ~0~ & ~0~\\ &\cdot&\\ &\cdot&\\ &\cdot&\\ ~0~ & ~0~ & ~0~\\ ~0~ & ~0~ & ~0~\\ &\cdot&\\ &\cdot&\\ &\cdot&\\ ~0~ & ~0~ & ~0~\\ ~1~ & ~x^1~ & ~-1~\\ &\cdot&\\ &\cdot&\\ &\cdot&\\ ~1~ & ~x^K~ & ~-1~\\ ~1~ & ~x^1~ & ~0~\\ &\cdot&\\ &\cdot&\\ &\cdot&\\ ~1~ & ~x^K~ & ~0~ \end{array} \right].(StopMathJax)

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).

negLogLik.m117 d00 = rho*capPi*diag(pExit(:,1));
negLogLik.m118 d01 = rho*capPi*diag(pExit(:,2));
negLogLik.m119 d10 = rho*capPi-d00;
negLogLik.m120 d11 = rho*capPi-d01;
negLogLik.m121 dPsi_dUbar = [[d00;d00;zeros(2*nSuppX,nSuppX)] [zeros(2*nSuppX,nSuppX);d01;d01]
negLogLik.m122 [d10;d10;zeros(2*nSuppX,nSuppX)] [zeros(2*nSuppX,nSuppX);d11;d11]];

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).

negLogLik.m127 dPsi_dTheta = [[zeros(2*nSuppX,1);ones(2*nSuppX,1)] [zeros(2*nSuppX,1);supportX;supportX] [zeros(2*nSuppX,1);-ones(nSuppX,1);zeros(nSuppX,1)]];

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

negLogLik.m131 dUbar_dTheta = (eye(4*nSuppX)-dPsi_dUbar)\dPsi_dTheta;
negLogLik.m132 dDeltaU_dTheta = dUbar_dTheta(2*nSuppX+1:4*nSuppX,:)-dUbar_dTheta(1:2*nSuppX,:);

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.

negLogLik.m136 nTheta = size(dUbar_dTheta,2);
negLogLik.m137 negFirmScores = repmat((1-2*choices).*(1-p),[1 1 nTheta]);
negLogLik.m138 for i=1:nTheta
negLogLik.m139 negFirmScores(:,:,i) = negFirmScores(:,:,i).*dDeltaU_dTheta(iX+nSuppX*laggedChoices+2*(i-1)*nSuppX);
negLogLik.m140 end
negLogLik.m141 negFirmScores=squeeze(sum(negFirmScores,1));
negLogLik.m142 negScore = sum(negFirmScores)';
negLogLik.m143end

Information Matrix

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

(0)

(StartMathJax) \sum_{n=1}^N\frac{\partial\log\left[\prod_{t=1}^Tp(a_{tn}|x_{tn},a_{(t-1)n})\right]}{\partial\theta}\cdot\frac{\partial\log\left[\prod_{t=1}^Tp(a_{tn}|x_{tn},a_{(t-1)n})\right]}{\partial\theta'}, (StopMathJax)

and returns it in informationMatrix.

negLogLik.m152if nargout==3
negLogLik.m153 informationMatrix = zeros(nTheta,nTheta);
negLogLik.m154 for n=1:size(negFirmScores,1)
negLogLik.m155 informationMatrix = informationMatrix + negFirmScores(n,:)'*negFirmScores(n,:);
negLogLik.m156 end
negLogLik.m157end

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).

5The Script that Puts It All Together

The script in dynamicDiscreteChoice.m simulates data and computed maximum likelihood estimates using the nested fixed point (NFXP) method of Rust (1987) and Rust (1994). It takes (StartMathJax){\cal X},\delta_0,\rho(StopMathJax) to be known, either takes (StartMathJax)\Pi(StopMathJax) to be known or estimates it in a first stage, and focuses on maximum partial likelihood estimation of the remaining parameters; (StartMathJax)\beta_0,\beta_1,\delta_1(StopMathJax); from conditional choice probabilities.

51Simulating Data

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

dynamicDiscreteChoice.m76nPeriods = 100
dynamicDiscreteChoice.m77nFirms = 1000

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.

dynamicDiscreteChoice.m81tolFixedPoint = 1e-10

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).

dynamicDiscreteChoice.m93nSuppX = 5;
dynamicDiscreteChoice.m94supportX = (1:nSuppX)'
dynamicDiscreteChoice.m95capPi = 1./(1+abs(ones(nSuppX,1)*(1:nSuppX)-(1:nSuppX)'*ones(1,nSuppX)));
dynamicDiscreteChoice.m96capPi = capPi./(sum(capPi')'*ones(1,nSuppX))
dynamicDiscreteChoice.m97beta = [-0.1*nSuppX;0.2]
dynamicDiscreteChoice.m98delta = [0;1]
dynamicDiscreteChoice.m99rho = 0.95

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).

dynamicDiscreteChoice.m103[u0,u1] = flowpayoffs(supportX,beta,delta);
dynamicDiscreteChoice.m104[capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,@bellman,[],[]);
dynamicDiscreteChoice.m105deltaU = capU1-capU0;

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).

dynamicDiscreteChoice.m109[choices,iX] = simulateData(deltaU,capPi,nPeriods,nFirms);

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.

dynamicDiscreteChoice.m115objectiveFunction = @(parameters)negLogLik(choices,iX,supportX,capPi,parameters(1:2),[delta(1);parameters(3)],
dynamicDiscreteChoice.m116 rho,@flowpayoffs,@bellman,@fixedPoint,tolFixedPoint)

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).

dynamicDiscreteChoice.m120startvalues = [-1;-0.1;0.5];

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.2

dynamicDiscreteChoice.m124lowerBounds = -Inf*ones(size(startvalues));
dynamicDiscreteChoice.m125lowerBounds(3) = 0;

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 installed3).

dynamicDiscreteChoice.m129OptimizerOptions = optimset('Display','iter','Algorithm','interior-point','AlwaysHonorConstraints','bounds',
dynamicDiscreteChoice.m130 'GradObj','on','TolFun',1E-6,'TolX',1E-10,'DerivativeCheck','off','TypicalX',[beta;delta(2)]);
dynamicDiscreteChoice.m131[maxLikEstimates,~,exitflag] = fmincon(objectiveFunction,startvalues,[],[],[],[],lowerBounds,[],[],OptimizerOptions)

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.

dynamicDiscreteChoice.m135[~,~,informationMatrix] = objectiveFunction(maxLikEstimates);
dynamicDiscreteChoice.m136standardErrors = diag(sqrt(inv(informationMatrix)));

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).

dynamicDiscreteChoice.m140disp('Summary of Results');
dynamicDiscreteChoice.m141disp('--------------------------------------------');
dynamicDiscreteChoice.m142disp(' true start estim ste.');
dynamicDiscreteChoice.m143disp([[beta;delta(2)] startvalues maxLikEstimates standardErrors]);

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.

dynamicDiscreteChoice.m148piHat = estimatePi(iX,nSuppX)

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.

6Student 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.


1The code has been tested with Matlab 2014a and later, with its Optimization Toolbox, on iMacs with OS X 10.9 and later. A free trial version of Knitro is available from Artelys.
2Note that fmincon, but also its alternatives discussed below, allow the user to specify bounds on parameters; if another function is used that does not allow for bounds on the parameters, you can use an alternative parameterization to ensure that parameters only take values in some admissible set (for example, you can specify (StartMathJax)\delta_1=\exp(\delta_1^*)(StopMathJax) for (StartMathJax)\delta_1^*\in\mathbb{R}(StopMathJax) to ensure that (StartMathJax)\delta_1 \gt 0(StopMathJax)). Minimizers like fmincon also allow you to impose more elaborate constraints on the parameters; you will need this option when implementing the MPEC alternative to NFXP of Judd and Su (2012) (see Section 6).
3fmincon requires Matlab's Optimization Toolbox, knitromatlab is included in Knitro 9.0, knitrolink uses both, and ktrlink can be used if the Optimization Toolbox is installed with an earlier version of Knitro.

1Contractions and Dynamic Programming

See e.g. Stokey, Lucas, and Prescott (1989) and Judd (1998) for a rigorous and comprehensive treatment of the topics in this appendix and their economic applications.

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

(0)

(StartMathJax)U_{n+1}=U_n-\left[I-\Psi'_{U_n}\right]^{-1}\left[U_n-\Psi(U_n)\right],(StopMathJax)

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).

2Miscellaneous Utilities

This section documents some miscellaneous utilities that are called by the main functions. At this point, it only includes randomDiscrete.

The function randomDiscrete returns a random draw from the distribution of (StartMathJax)(Y_1,\ldots,Y_n)(StopMathJax); with (StartMathJax)Y_1,\ldots,Y_n(StopMathJax) independently discretely distributed with (not necessarily identical) distributions on (StartMathJax)\{1,2,\ldots,k\}(StopMathJax) for some (StartMathJax)2\leq k \lt \infty(StopMathJax).

randomDiscrete.m7function y = randomDiscrete(p)

It requires one input argument:

{p} a (StartMathJax)k\times n(StopMathJax) matrix with (StartMathJax)(k,n)(StopMathJax)th element (StartMathJax)\Pr(Y_n=k)(StopMathJax).

It returns

{y} a (StartMathJax)1\times n(StopMathJax) vector with a random draw (StartMathJax)(y_1,\ldots,y_N)(StopMathJax) from (StartMathJax)(Y_1,\ldots,Y_N)(StopMathJax).

The function randomDiscrete first stores the number (StartMathJax)k(StopMathJax) of support points in a scalar nSupp and the number (StartMathJax)n(StopMathJax) of random variables in a scalar nVar.

randomDiscrete.m19nSupp = size(p,1);
randomDiscrete.m20nVar = size(p,2);

Then, it creates a (StartMathJax)(k-1)\times n(StopMathJax) matrix uniformDraws with (StartMathJax)k-1(StopMathJax) identical rows, containing (StartMathJax)n(StopMathJax) independent draws from a standard uniform distribution, and computes the (StartMathJax)k\times n(StopMathJax) matrix cumulativeP with (StartMathJax)(k,n)(StopMathJax)th element (StartMathJax)\Pr(Y_n\leq k)(StopMathJax).

randomDiscrete.m24uniformDraws = ones(nSupp-1,1)*random('unif',zeros(1,nVar),ones(1,nVar));
randomDiscrete.m25cumulativeP = cumsum(p);

Finally, for each (StartMathJax)n(StopMathJax), it sets (StartMathJax)y_n(StopMathJax) equal to 1 plus the number of elements of (StartMathJax)\{\Pr(Y_n\leq 1), \ldots, \Pr(Y_n\leq k-1)(StopMathJax)} that are weakly smaller than the uniform random draw in the (StartMathJax)n(StopMathJax)th column of uniformDraws.

randomDiscrete.m29y = sum([ones(1,nVar);cumulativeP(1:nSupp-1,:)<=uniformDraws]);

dynamicDiscreteChoice.m76nPeriods = 100
dynamicDiscreteChoice.m77nFirms = 1000
dynamicDiscreteChoice.m81tolFixedPoint = 1e-10
dynamicDiscreteChoice.m93nSuppX = 5;
dynamicDiscreteChoice.m94supportX = (1:nSuppX)'
dynamicDiscreteChoice.m95capPi = 1./(1+abs(ones(nSuppX,1)*(1:nSuppX)-(1:nSuppX)'*ones(1,nSuppX)));
dynamicDiscreteChoice.m96capPi = capPi./(sum(capPi')'*ones(1,nSuppX))
dynamicDiscreteChoice.m97beta = [-0.1*nSuppX;0.2]
dynamicDiscreteChoice.m98delta = [0;1]
dynamicDiscreteChoice.m99rho = 0.95
dynamicDiscreteChoice.m103[u0,u1] = flowpayoffs(supportX,beta,delta);
dynamicDiscreteChoice.m104[capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,@bellman,[],[]);
dynamicDiscreteChoice.m105deltaU = capU1-capU0;
dynamicDiscreteChoice.m109[choices,iX] = simulateData(deltaU,capPi,nPeriods,nFirms);
dynamicDiscreteChoice.m115objectiveFunction = @(parameters)negLogLik(choices,iX,supportX,capPi,parameters(1:2),[delta(1);parameters(3)],
dynamicDiscreteChoice.m116 rho,@flowpayoffs,@bellman,@fixedPoint,tolFixedPoint)
dynamicDiscreteChoice.m120startvalues = [-1;-0.1;0.5];
dynamicDiscreteChoice.m124lowerBounds = -Inf*ones(size(startvalues));
dynamicDiscreteChoice.m125lowerBounds(3) = 0;
dynamicDiscreteChoice.m129OptimizerOptions = optimset('Display','iter','Algorithm','interior-point','AlwaysHonorConstraints','bounds',
dynamicDiscreteChoice.m130 'GradObj','on','TolFun',1E-6,'TolX',1E-10,'DerivativeCheck','off','TypicalX',[beta;delta(2)]);
dynamicDiscreteChoice.m131[maxLikEstimates,~,exitflag] = fmincon(objectiveFunction,startvalues,[],[],[],[],lowerBounds,[],[],OptimizerOptions)
dynamicDiscreteChoice.m135[~,~,informationMatrix] = objectiveFunction(maxLikEstimates);
dynamicDiscreteChoice.m136standardErrors = diag(sqrt(inv(informationMatrix)));
dynamicDiscreteChoice.m140disp('Summary of Results');
dynamicDiscreteChoice.m141disp('--------------------------------------------');
dynamicDiscreteChoice.m142disp(' true start estim ste.');
dynamicDiscreteChoice.m143disp([[beta;delta(2)] startvalues maxLikEstimates standardErrors]);
dynamicDiscreteChoice.m148piHat = estimatePi(iX,nSuppX)
flowpayoffs.m7function [u0,u1] = flowpayoffs(supportX,beta,delta)
flowpayoffs.m24nSuppX = size(supportX,1);
flowpayoffs.m55u0 = [zeros(nSuppX,1) -delta(1)*ones(nSuppX,1)];
flowpayoffs.m56u1 = [ones(nSuppX,1) supportX]*beta*[1 1]-delta(2)*ones(nSuppX,1)*[1 0];
bellman.m15function [capU0,capU1] = bellman(capU0,capU1,u0,u1,capPi,rho)
bellman.m34r0 = log(exp(capU0(:,1))+exp(capU1(:,1)));
bellman.m35r1 = log(exp(capU0(:,2))+exp(capU1(:,2)));
bellman.m39capU0 = u0 + rho*capPi*r0*[1 1];
bellman.m40capU1 = u1 + rho*capPi*r1*[1 1];
fixedPoint.m7function [capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,bellman,capU0,capU1)
fixedPoint.m28nSuppX = size(capPi,1);
fixedPoint.m32if isempty(capU0)
fixedPoint.m33 capU0 = zeros(nSuppX,2);
fixedPoint.m34end
fixedPoint.m35if isempty(capU1)
fixedPoint.m36 capU1 = zeros(nSuppX,2);
fixedPoint.m37end
fixedPoint.m41inU0 = capU0+2*tolFixedPoint;
fixedPoint.m42inU1 = capU1+2*tolFixedPoint;
fixedPoint.m43while (max(max(abs(inU0-capU0)))>tolFixedPoint) || (max(max(abs(inU1-capU1)))>tolFixedPoint);
fixedPoint.m44 inU0 = capU0;
fixedPoint.m45 inU1 = capU1;
fixedPoint.m46 [capU0,capU1] = bellman(inU0,inU1,u0,u1,capPi,rho);
fixedPoint.m47end
simulateData.m7function [choices,iX] = simulateData(deltaU,capPi,nPeriods,nFirms)
simulateData.m23nSuppX = size(capPi,1);
simulateData.m39oneMinPi = eye(nSuppX)-capPi';
simulateData.m40pInf = [oneMinPi(1:nSuppX-1,:);ones(1,nSuppX)]\[zeros(nSuppX-1,1);1];
simulateData.m44iX = randomDiscrete(pInf*ones(1,nFirms));
simulateData.m48deltaEpsilon = random('ev',zeros(1,nFirms),ones(1,nFirms))-random('ev',zeros(1,nFirms),ones(1,nFirms));
simulateData.m49choices = deltaU(iX,1)' > deltaEpsilon;
simulateData.m53for t = 2:nPeriods
simulateData.m54 iX = [iX;randomDiscrete(capPi(iX(end,:),:)')];
simulateData.m55 deltaEpsilon = random('ev',zeros(1,nFirms),ones(1,nFirms))-random('ev',zeros(1,nFirms),ones(1,nFirms));
simulateData.m56 choices = [choices;(deltaU(iX(end,:)+nSuppX*choices(end,:)) > deltaEpsilon)];
simulateData.m57end
estimatePi.m7function piHat = estimatePi(iX,nSuppX)
estimatePi.m20nPeriods=size(iX,1);
estimatePi.m24for i=1:nSuppX
estimatePi.m25 for j=1:nSuppX
estimatePi.m26 piHat(i,j)=sum(sum((iX(2:nPeriods,:)==j)&(iX(1:nPeriods-1,:)==i)))/sum(sum((iX(1:nPeriods-1,:)==i)));
estimatePi.m27 end
estimatePi.m28end
negLogLik.m7function [nll,negScore,informationMatrix] =
negLogLik.m8 negLogLik(choices,iX,supportX,capPi,beta,delta,rho,flowpayoffs,bellman,fixedPoint,tolFixedPoint)
negLogLik.m35nSuppX = size(supportX,1);
negLogLik.m39[u0,u1] = flowpayoffs(supportX,beta,delta);
negLogLik.m40[capU0,capU1] = fixedPoint(u0,u1,capPi,rho,tolFixedPoint,bellman,[],[]);
negLogLik.m41deltaU = capU1-capU0;
negLogLik.m42pExit = 1./(1+exp(deltaU));
negLogLik.m49laggedChoices = [zeros(1,size(choices,2));choices(1:end-1,:)];
negLogLik.m50p = choices + (1-2*choices).*pExit(iX+nSuppX*laggedChoices);
negLogLik.m51nll = -sum(sum(log(p)));
negLogLik.m57if nargout>=2
negLogLik.m117 d00 = rho*capPi*diag(pExit(:,1));
negLogLik.m118 d01 = rho*capPi*diag(pExit(:,2));
negLogLik.m119 d10 = rho*capPi-d00;
negLogLik.m120 d11 = rho*capPi-d01;
negLogLik.m121 dPsi_dUbar = [[d00;d00;zeros(2*nSuppX,nSuppX)] [zeros(2*nSuppX,nSuppX);d01;d01]
negLogLik.m122 [d10;d10;zeros(2*nSuppX,nSuppX)] [zeros(2*nSuppX,nSuppX);d11;d11]];
negLogLik.m127 dPsi_dTheta = [[zeros(2*nSuppX,1);ones(2*nSuppX,1)] [zeros(2*nSuppX,1);supportX;supportX] [zeros(2*nSuppX,1);-ones(nSuppX,1);zeros(nSuppX,1)]];
negLogLik.m131 dUbar_dTheta = (eye(4*nSuppX)-dPsi_dUbar)\dPsi_dTheta;
negLogLik.m132 dDeltaU_dTheta = dUbar_dTheta(2*nSuppX+1:4*nSuppX,:)-dUbar_dTheta(1:2*nSuppX,:);
negLogLik.m136 nTheta = size(dUbar_dTheta,2);
negLogLik.m137 negFirmScores = repmat((1-2*choices).*(1-p),[1 1 nTheta]);
negLogLik.m138 for i=1:nTheta
negLogLik.m139 negFirmScores(:,:,i) = negFirmScores(:,:,i).*dDeltaU_dTheta(iX+nSuppX*laggedChoices+2*(i-1)*nSuppX);
negLogLik.m140 end
negLogLik.m141 negFirmScores=squeeze(sum(negFirmScores,1));
negLogLik.m142 negScore = sum(negFirmScores)';
negLogLik.m143end
negLogLik.m152if nargout==3
negLogLik.m153 informationMatrix = zeros(nTheta,nTheta);
negLogLik.m154 for n=1:size(negFirmScores,1)
negLogLik.m155 informationMatrix = informationMatrix + negFirmScores(n,:)'*negFirmScores(n,:);
negLogLik.m156 end
negLogLik.m157end
randomDiscrete.m7function y = randomDiscrete(p)
randomDiscrete.m19nSupp = size(p,1);
randomDiscrete.m20nVar = size(p,2);
randomDiscrete.m24uniformDraws = ones(nSupp-1,1)*random('unif',zeros(1,nVar),ones(1,nVar));
randomDiscrete.m25cumulativeP = cumsum(p);
randomDiscrete.m29y = sum([ones(1,nVar);cumulativeP(1:nSupp-1,:)<=uniformDraws]);