Viewing and Using this File
This file (DynamicDiscreteChoice.m.html) 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 in 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.
This file was generated by applying the Comments++ tool of Campbell and Peters (2010) to the Matlab script DynamicDiscreteChoice.m and the functions it calls. This script and these functions are all written as “literate programs” that embed their own documentation in a formatting language similar to LaTeX.
Figure 1: An Example of Displayed Code from this Document
Below, working code is displayed in areas formatted with “bars” similar to those on “old-style” computer paper for dot-matrix printers. Figure 1 shows a screen shot of a few lines of code from this document. The source file name is to the left, and the line number from the source file preceeds each line. In this example, Line 54 is highlighted blue because the mouse was hovering over it when screen shot was taken.
Figure 2: Switching Styles in Firefox
You can change this document to show only the working code. Figure 2 shows how to do this in Firefox.1 The default value for Page Style is Display All. Changing this to Display Only Code hides the documentation.
Printing this File
So that the displayed line numbers match with the actual line numbers, your browser preserves all line breaks in the displayed source code. That is, a long line cannot break into multiple lines after you resize your browser window. The same principle applies to printing. If a line of code is too long to be displayed on the page, its right end will be cut off. For this reason, this document is best printed in landscape mode.2
Software Requirements
Running the code documented in this file requires Matlab version 2008a or up, with its Optimization Toolbox and, ideally, the Knitro libraries. You can substitute a standard optimization function from the Optimization Toolbox, like fmincon, for Knitro's ktrlink if you do not have access to the Knitro libraries. However, note that a basic student version of these libraries is available from Ziena Optimization for free.
Road Map
The remainder of this document proceeds as follows. The next section covers two functions that define the decision problem, flowpayoffs and Psi. 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 2B discusses FixedPointPsi, 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 defined by Psi. Section 3C documents how SimulateDP can be used to simulate data for the decision problem set up by flowpayoffs and Psi. Section 4D documents functions for estimating the transition matrix 
and for computing the log partial likelihood for conditional choice data. Section 5E 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 6F contains a range of student exercises. Appendix 1A provides mathematical background by discussing some theory of contractions. Appendix 2B documents procedures that are called by the main routines, but that are of limited substantial interest, such as RandomDiscr.
1AA 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 
, 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 
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, 
). They also depend on externally determined (choice-independent) observed scalar state variables 
, as well as unobserved state variables 
and 
that are revealed to the firm before it makes its choice in period 
. Specifically, in period 
, its flow profits from choice 0 are
(0)
and its flow profits from choice 1 are
(0)
Note that 
and 
are exit and entry costs that the firm only pays if it changes state. Gross of these costs, an active firm makes profits 
and an inactive firm has profits 
.
The state variable 
has finite support 
. From its random initial value 
, it follows a first-order Markov chain with 
transition probability matrix 
, with typical element 
, independently of the firm's choices. The profits shocks 
are independent of 
, across time 
, and choices 
, with type I extreme value distributions centered around 0. Note that the firm controls the evolution of the state 
only through its choice of 
, as 
and 
are “exogenous” state variables of which the evolution is externally specified.
The firm has rational expectations about future states and choices. It chooses the action 
that maximizes its expected flow of profits, discounted at a factor 
.
Two functions, flowpayoffs and Psi, 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 Psi.m.
1A.1Flow Payoffs
The function flowpayoffs computes the mean (over 
) flow payoffs, 
and 
, for each profit and past choice pair 
.
| flowpayoffs.m | 7 | function [u0,u1] = flowpayoffs(calX,beta,delta) |
It requires the following input arguments:
- calX
- a
vector with the support points of the profit state 
(the elements of 
, consistently ordered with the Markov transition matrix 
);
- beta
- a
vector that contains the intercept ( 
) and profit state slope ( 
) of the net payoffs to choice 
; and
- delta
- a
vector that contains the firm's exit ( 
) and entry ( 
) costs.
It returns
- u0
- a
matrix of which the 
th entry is 
and
- u1
- a
matrix of which the 
th entry is 
.
That is, the rows correspond to the support points of 
, and the columns to the choice in the previous period, 
.
The function flowpayoffs first stores the number 
of elements of calX in a scalar K.
| flowpayoffs.m | 24 | K = size(calX,1); |
Then, it constructs a 
matrix u0 with the value of
(0)
and a 
matrix u1 with the value of
(0)
| flowpayoffs.m | 55 | u0 = [zeros(K,1) -delta(1)*ones(K,1)]; |
| flowpayoffs.m | 56 | u1 = [ones(K,1) calX]*beta*[1 1]-delta(2)*ones(K,1)*[1 0]; |
You can change the specification of the flow profits by editing these two lines of code.
1A.2Bellman Operator
The expected discounted profits net of 
immediately following choice 
at time 
, 
, satisfy a recursive system 
and 
or, with 
and 
, simply 
(see e.g. Rust (1994)). Here, 
is a Bellman-like operator that embodies the model's dynamic specification and that depends on the flow payoffs 
and 
, the Markov transition matrix 
, and the discount factor 
. Its elements are the operators 
and 
, with 
implicitly defined by the right hand side of
(1)
where 
is McFadden (1981)'s social surplus for the binary choice problem with utilities 
and 
. The first term in the right-hand side of (1) equals period 
's flow profits following choice 
, net of 
. The second term equals the expected discounted profits from continuing into period 
with choice 
(note that, given the choice 
, these continuation payoffs do not depend on the past choice 
). The (mean-zero) extreme value assumptions on 
and 
implies that
(2)
The function Psi applies the operator 
once to input values of 
and returns 
.
| Psi.m | 15 | function [capU0,capU1] = Psi(capU0,capU1,u0,u1,Pi,rho) |
It requires the following input arguments:
- capU0
- a
matrix of which the 
th entry is 
;
- capU1
- a
matrix of which the 
th entry is 
;
- u0
- a
matrix of which the 
th entry is 
;
- u1
- a
matrix of which the 
th entry is 
;
- Pi
- the
Markov transition matrix 
for 
, with typical element 
; and
- rho
- a scalar with the value of the discount factor
.
It returns
- capU0
- a
matrix of which the 
th entry is 
;
- capU1
- a
matrix of which the 
th entry is 
;
To this end, Psi first computes the surpluses 
and 
in (2) for all 
and stacks these in 
vectors R0 and R1.
| Psi.m | 33 | R0 = log(exp(capU0(:,1))+exp(capU1(:,1))); |
| Psi.m | 34 | R1 = log(exp(capU0(:,2))+exp(capU1(:,2))); |
Then, it applies (1) to compute new values of capU0 and capU1.
| Psi.m | 38 | capU0 = u0 + rho*Pi*R0*[1 1]; |
| Psi.m | 39 | capU1 = u1 + rho*Pi*R1*[1 1]; |
Here, the conditional expectation over 
in (1) is taken by premultiplying the vectors R0 and R1 by the Markov transition matrix Pi. 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 
, and therefore the function Psi, through the specification of the surpluses 
and 
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).
2BSolving the Decision Problem
The function FixedPointPsi computes the fixed point 
of the Bellman-like operator 
, using the method of successive approximations. It is easy to show that 
is a contraction, so that the method of successive approximations converges linearly and globally to a point within a positive
maximum absolute distance tFP (see Appendix 1A.3).
| FixedPointPsi.m | 8 | function [capU0,capU1] = FixedPointPsi(u0,u1,Pi,rho,tFP,Psi,capU0,capU1) |
It requires the following input arguments:
- u0
- a
matrix of which the 
th entry is 
;
- u1
- a
matrix of which the 
th entry is 
;
- Pi
- the
Markov transition matrix 
for 
, with typical element 
;
- rho
- a scalar with the value of the discount factor
;
- tFP
- a scalar tolerance level that is used to determine convergence of the successive approximations;
- Psi
- the handle of the function [U0,U1]=Psi(Uin0,Uin1,u0,u1,Pi,rho) that iterates once on
;
- capU0
- a
matrix of which the 
th entry is a starting value for 
(optional; set to [] to select default starting value);
- capU1
- a
matrix of which the 
th entry is a starting value for 
(optional; set to [] to select default starting value);
It returns
- capU0
- a
matrix of which the 
th entry is 
; and
- capU1
- a
matrix of which the 
th entry is 
;
The function FixedPointPsi first stores the number 
of elements of calX in a scalar K.
| FixedPointPsi.m | 28 | K = size(Pi,1); |
The starting values for 
and 
are set to 
if the input arguments capU0 and capU1 are empty.
| FixedPointPsi.m | 32 | if isempty(capU0) |
| FixedPointPsi.m | 33 | capU0 = zeros(K,2); |
| FixedPointPsi.m | 34 | end |
| FixedPointPsi.m | 35 | if isempty(capU1) |
| FixedPointPsi.m | 36 | capU1 = zeros(K,2); |
| FixedPointPsi.m | 37 | end |
The 
matrices Uin0 and Uin1 store the value of 
that is fed into the operator 
, for comparison with the value of 
that is subsequently stored in capU0 and capU1. They are initialized to deviate from capU0 and capU1 by more than tFP, so that the while statement allows at least one iteration of 
, and stops as soon as 
no longer exceeds the tolerance level in tFP.
| FixedPointPsi.m | 41 | Uin0 = capU0+2*tFP; |
| FixedPointPsi.m | 42 | Uin1 = capU1+2*tFP; |
| FixedPointPsi.m | 43 | while (max(max(abs(Uin0-capU0)))>tFP) || (max(max(abs(Uin1-capU1)))>tFP); |
| FixedPointPsi.m | 44 | Uin0 = capU0; |
| FixedPointPsi.m | 45 | Uin1 = capU1; |
| FixedPointPsi.m | 46 | [capU0,capU1] = Psi(Uin0,Uin1,u0,u1,Pi,rho); |
| FixedPointPsi.m | 47 | end |
You can replace FixedPointPsi by another function if you want to use alternative methods, such as Newton methods, for computing the fixed point 
, or if you want to work with finite horizon problems.
3CData Simulation
Suppose that we have computed 
for all 
, using e.g. flowpayoffs and FixedPointPsi, and that we have specified the Markov transition matrix 
. Then, the function SimulateDP can be used to simulate 
independent histories 
, taking 
as given and drawing 
from the stationary distribution of 
.
| SimulateDP.m | 7 | function [A,iX] = SimulateDP(DeltaU,Pi,T,N) |
The function SimulateDP requires the following input arguments:
- DeltaU
- a
matrix of which the 
th entry is 
;
- Pi
- the
Markov transition matrix 
for 
, with typical element 
;
- T
- the scalar number
of time periods to simulate data for; and
- N
- the scalar number
of firms to simulate data for.
It returns
- A
- a
matrix with simulated choices, with each each column containing an independent simulation of 
; and
- iX
- a
matrix with simulated states, with each each column containing the indices (in 
) of an independent simulation of 
. For example, if 
for some firm, then 3, not 
, is stored in iX.
The function SimulateDP first stores the number 
of elements of calX in a scalar K.
| SimulateDP.m | 23 | K = size(Pi,1); |
Next, it assumes that 
is ergodic and initializes the simulation by drawing 
, 
independent times, from the stationary distribution of 
. To this end, it first computes the 
vector 
with the stationary probabilities 
using
(0)
and stores it in Pinf.
| SimulateDP.m | 39 | oneminPi = eye(K)-Pi'; |
| SimulateDP.m | 40 | Pinf = inv([oneminPi(1:K-1,:);ones(1,K)])*[zeros(K-1,1);1]; |
Then, it uses the auxiliary function RandomDiscr (see Appendix 2B) and the values stored in Pinf to simulate a 
vector of values of 
from the stationary distribution 
and stores their indices in iX.
| SimulateDP.m | 44 | iX = RandomDiscr(Pinf*ones(1,N)); |
Using these 
simulated values of 
, and 
simulated values of 
which are stored in Depsilon, it simulates 
values of the first choice by using that 
if 
and 
otherwise. These are stored in the 
vector A.
| SimulateDP.m | 48 | Depsilon = random('ev',zeros(1,N),ones(1,N))-random('ev',zeros(1,N),ones(1,N)); |
| SimulateDP.m | 49 | A = DeltaU(iX,1)' > Depsilon; |
Finally, 
values of 
are simulated, using the transition matrix 
and RandomDiscr, and their indices added as a row to the bottom of the 
matrix iX; and 
values of 
are simulated, using that 
if 
and 
otherwise, and stored as a row at the bottom of the 
matrix A; recursively for 
.
| SimulateDP.m | 53 | for t = 2:T |
| SimulateDP.m | 54 | iX = [iX;RandomDiscr(Pi(iX(end,:),:)')]; |
| SimulateDP.m | 55 | Depsilon = random('ev',zeros(1,N),ones(1,N))-random('ev',zeros(1,N),ones(1,N)); |
| SimulateDP.m | 56 | A = [A;(DeltaU(iX(end,:)+K*A(end,:)) > Depsilon)]; |
| SimulateDP.m | 57 | end |
4DEstimation
Suppose that we have a random sample 
from the population of state and choice histories 
. Following Rust (1994), we provide functions for implementing a two-stage procedure in which, first, the state transition matrix 
is estimated directly from observed state transitions and, second, the remaining parameters are estimated by maximizing the partial likelihood for the observed choices.
4D.1First Stage: State Transitions
The function estimatePi computes a frequency estimate of the Markov transition matrix 
from state transition data.
| estimatePi.m | 7 | function Pihat = estimatePi(iX,K) |
It requires the following input arguments:
- iX
- a
matrix with indices of observed states 
in 
(for example, if 
, then the first element of iX is 3, not 
);
- K
- the scalar number
of support points of the profit state 
(the number of elements of 
)
and returns
- Pihat
- an estimate
of the 
Markov transition matrix 
for 
, with typical element 
equal to the sample frequency of transitions to state 
among all transitions from state 
.
The function estimatePi first stores the number of time periods 
in a scalar T.
| estimatePi.m | 20 | T=size(iX,1); |
Then, for each pair 
, it estimates the probability 
by the appropriate sample frequency, the number of transitions from 
to 
divided by the total number of transitions from 
in the data iX.
| estimatePi.m | 24 | for i=1:K |
| estimatePi.m | 25 | for j=1:K |
| estimatePi.m | 26 | Pihat(i,j)=sum(sum((iX(2:T,:)==j)&(iX(1:T-1,:)==i)))/sum(sum((iX(1:T-1,:)==i))); |
| estimatePi.m | 27 | end |
| estimatePi.m | 28 | end |
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 smooths across observations.
4D.2Second Stage: Choices
The function NegLogLikCCP computes minus the log partial likelihood for the conditional choice part of the data. Optionally, it also returns the corresponding score vector and an estimate of the information matrix for the parameter (sub)vector 
(the scores are specific to the estimation example in Section 5E's script and should be adapted for inference on other parameters).
| NegLogLikCCP.m | 7 | function [nll grad gradgrad] = NegLogLikCCP(A,iX,calX,Pi,beta,delta,rho,flowpayoffs,Psi,FixedPointPsi,tFP) |
The function NegLogLikCCP requires the following input arguments:
- A
- a
matrix with choice observations 
;
- iX
- a
matrix with indices of observed states 
in 
(for example, if 
, then the first element of iX is 3, not 
);
- calX
- a
vector with the support points of the profit state 
(the elements of 
, consistently ordered with the Markov transition matrix 
);
- Pi
- the (possibly estimated)
Markov transition matrix 
for 
, with typical element 
;
- beta
- a
vector that contains the intercept ( 
) and profit state slope ( 
) of the flow payoffs to choice 
;
- delta
- a
vector that contains the firm's exit ( 
) and entry ( 
) costs.
- rho
- a scalar with the value of the discount factor
;
- flowpayoffs
- a handle of a function [u0,u1]=flowpayoffs(calX,beta,delta) that computes the mean flow payoffs
and 
;
- Psi
- a handle of a function [U0,U1]=Psi(Uin0,Uin1,u0,u1,Pi,rho) that iterates once on
;
- FixedPointPsi
- a handle of a function [U0,U1]=FixedPointPsi(u0,u1,Pi,rho,tFP,Psi,U0,U1) that computes the fixed point
of 
; and
- tFP
- a scalar tolerance level that is used to determine convergence of the successive approximations of the fixed point
of 
.
It returns
- nll
- a scalar with minus the log partial likelihood for the conditional choices
and optionally
- grad
- a
vector with minus the partial likelihood score for 
and
- gradgrad
- a
matrix with the sum of the 
outer products of the individual contributions to the score for 
.
The function NegLogLikCCP first stores the number 
of elements of calX in a scalar K.
| NegLogLikCCP.m | 34 | K = size(calX,1); |
Next, it computes the flow payoffs 
(u0) and 
(u1), the choice-specific net expected discounted values 
(capU0) and 
(capU1), their contrast 
(DeltaU), and the implied probabilities 
of not serving the market (pexit) for the inputted parameter values. Note that this implements the NFXP procedure's inner loop.
| NegLogLikCCP.m | 38 | [u0,u1] = flowpayoffs(calX,beta,delta); |
| NegLogLikCCP.m | 39 | [capU0,capU1] = FixedPointPsi(u0,u1,Pi,rho,tFP,Psi,[],[]); |
| NegLogLikCCP.m | 40 | DeltaU = capU1-capU0; |
| NegLogLikCCP.m | 41 | pexit = 1./(1+exp(DeltaU)); |
Log Partial Likelihood
The contribution to the likelihood of firm 
's choice in period 
is the conditional choice probability
(0)
with 
. The function NegLogLikCCP first computes these probabilities for each firm 
and period 
and stores them in a 
matrix p. Then, it returns minus the sum of their logs, the log partial likelihood for the conditional choices, in nll.
| NegLogLikCCP.m | 48 | Alag = [zeros(1,size(A,2));A(1:end-1,:)]; |
| NegLogLikCCP.m | 49 | p = A + (1-2*A).*pexit(iX+K*(Alag)); |
| NegLogLikCCP.m | 50 | nll = -sum(sum(log(p))); |
Score
If two or more output arguments are demanded from NegLogLikCCP, it computes and returns minus the partial likelihood
score for 
(the derivative of minus the log partial likelihood with respect to 
), in grad.
| NegLogLikCCP.m | 56 | if nargout>=2 |
Firm 
's contribution to the score equals
(0)
Its calculation requires that we compute 
. Recall that 
, and therewith 
, is only implicitly given by 
. Note that 
is defined on a set with 
points, so that 
can be represented by a 
vector and 
by a mapping from 
into 
. Specifically, 
can be represented by the 
vector 
that lists the values of 
and 
on their domain,
(0)
and that satisfies
(0)
with 
an appropriately rearranged version of 
(note that we made its dependence on 
explicit). With this alternative representation of 
and 
in place, 
follows from (see also Rust (1994), p.3110)
(0)
where 
is a 
unit matrix;
(0)
with
(0)
and 
is a 
matrix of zeros; and
(0)
The function NegLogLikCCP first computes 
(D00), 
(D01), 
(D10), 
(D11) and constructs 
(dPsi_dUbar).
| NegLogLikCCP.m | 115 | D00 = rho*Pi*diag(pexit(:,1)); |
| NegLogLikCCP.m | 116 | D01 = rho*Pi*diag(pexit(:,2)); |
| NegLogLikCCP.m | 117 | D10 = rho*Pi-D00; |
| NegLogLikCCP.m | 118 | D11 = rho*Pi-D01; |
| NegLogLikCCP.m | 119 | dPsi_dUbar = [[D00;D00;zeros(2*K,K)] [zeros(2*K,K);D01;D01] [D10;D10;zeros(2*K,K)] [zeros(2*K,K);D11;D11]]; |
It then computes 
(dPsi_dpar;
note that the following line is the only code in NegLogLikCCP that needs to be changed if other parameters that those in 
are estimated).
| NegLogLikCCP.m | 124 | dPsi_dpar = [[zeros(2*K,1);ones(2*K,1)] [zeros(2*K,1);calX;calX] [zeros(2*K,1);-ones(K,1);zeros(K,1)]]; |
Next, it computes 
(dUbar_dpar) and 
(dDeltaU_dpar).
| NegLogLikCCP.m | 128 | dUbar_dpar = inv(eye(4*K)-dPsi_dUbar)*dPsi_dpar; |
| NegLogLikCCP.m | 129 | dDeltaU_dpar = dUbar_dpar(2*K+1:4*K,:)-dUbar_dpar(1:2*K,:); |
Finally, it computes the 
vector 
for each 
, stacks these individual (minus) score contributions in the 
matrix indgrad, and sums them to compute minus the score vector, grad.
| NegLogLikCCP.m | 133 | npar = size(dUbar_dpar,2); |
| NegLogLikCCP.m | 134 | indgrad = repmat((1-2*A).*(1-p),[1 1 npar]); |
| NegLogLikCCP.m | 135 | for i=1:npar |
| NegLogLikCCP.m | 136 | indgrad(:,:,i) = indgrad(:,:,i).*dDeltaU_dpar(iX+K*Alag+2*(i-1)*K); |
| NegLogLikCCP.m | 137 | end |
| NegLogLikCCP.m | 138 | indgrad=squeeze(sum(indgrad,1)); |
| NegLogLikCCP.m | 139 | grad = sum(indgrad)'; |
| NegLogLikCCP.m | 140 | end |
Information Matrix
If three output arguments are demanded, NegLogLikCCP computes the sum of the 
outer products of the individual score contributions,
(0)
and returns it in gradgrad.
| NegLogLikCCP.m | 149 | if nargout==3 |
| NegLogLikCCP.m | 150 | gradgrad = zeros(npar,npar); |
| NegLogLikCCP.m | 151 | for n=1:size(indgrad,1) |
| NegLogLikCCP.m | 152 | gradgrad = gradgrad + indgrad(n,:)'*indgrad(n,:); |
| NegLogLikCCP.m | 153 | end |
| NegLogLikCCP.m | 154 | end |
When evaluated at an estimate of 
, this provides an estimate of the expected partial likelihood information matrix for 
. This estimate can in turn be used to estimate the variance-covariance matrix, and thus the standard errors, of the maximum partial likelihood estimator 
of 
.
5EThe 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 
to be known, either takes 
to be known or
estimates it in a first stage, and focuses on maximum partial likelihood estimation of the remaining parameters; 
; from conditional choice probabilities.
Simulating Data
First, we set the number of time periods (T) and firms (N) that we would like to have in our sample.
| DynamicDiscreteChoice.m | 93 | T = 100 |
| DynamicDiscreteChoice.m | 94 | N = 1000 |
We also set the tolerance tFP on the fixed point 
of 
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.m | 98 | tFP = 1e-10 |
Next, we specify the values of the model's parameters used in the simulation:
- K
- the scalar number
of elements of 
;
- calX
- the
vector 
with the support points of 
;
- Pi
- the
Markov transition matrix 
for 
, with typical element 
;
- beta
- the
vector 
with the parameters of the flow profit of active firms;
- delta
- the
vector of exit and entry costs 
;
- rho
- the scalar discount factor
;
| DynamicDiscreteChoice.m | 110 | K = 5; |
| DynamicDiscreteChoice.m | 111 | calX = (1:K)' |
| DynamicDiscreteChoice.m | 112 | Pi = 1./(1+abs(ones(K,1)*(1:K)-(1:K)'*ones(1,K))); |
| DynamicDiscreteChoice.m | 113 | Pi = Pi./(sum(Pi')'*ones(1,K)) |
| DynamicDiscreteChoice.m | 114 | beta = [-0.1*K;0.2] |
| DynamicDiscreteChoice.m | 115 | delta = [0;1] |
| DynamicDiscreteChoice.m | 116 | rho = 0.95 |
For these parameter values, we compute the flow payoffs 
(u0) and 
(u1), the choice-specific expected discounted values 
(U0) and 
(U1), and their contrast 
(DeltaU).
| DynamicDiscreteChoice.m | 120 | [u0,u1] = flowpayoffs(calX,beta,delta); |
| DynamicDiscreteChoice.m | 121 | [capU0,capU1] = FixedPointPsi(u0,u1,Pi,rho,tFP,@Psi,[],[]); |
| DynamicDiscreteChoice.m | 122 | DeltaU = capU1-capU0; |
With 
computed, and 
specified, we proceed to simulate a 
matrix of choices A and a 
matrix of states iX (recall from Section 3C that iX contains indices that point to elements of 
rather than those values themselves).
| DynamicDiscreteChoice.m | 126 | [A,iX] = SimulateDP(DeltaU,Pi,T,N); |
Nested Fixed Point Maximum Likelihood Estimation
First, set useKnitro to 1 if Knitro should be used for the outer (likelihood maximization) loop of the NFXP algorithm, and to 0 otherwise. Also, set the tolerances tc and tp for outer loop convergence (recall that the inner loop tolerance is already set in tFP).3
| DynamicDiscreteChoice.m | 132 | useKnitro = 0; |
| DynamicDiscreteChoice.m | 133 | tc = 1E-6 |
| DynamicDiscreteChoice.m | 134 | tp = 1E-6 |
Next, we specify a 
vector startvalues with starting values for the parameters to be estimated, 
.
| DynamicDiscreteChoice.m | 138 | startvalues = [-1;-0.1;0.5]; |
We continue with setting a lower bound (lb) of 0 on the third parameter, 
, and (nonbinding) lower bounds of 
on the other two parameters. There is no need to specify upper bounds.4
| DynamicDiscreteChoice.m | 142 | lb = -Inf*ones(size(startvalues)); |
| DynamicDiscreteChoice.m | 143 | lb(3) = 0; |
First, suppose that 
is known. We use NegLogLikCCP to specify an objective function objfun to be minimized in the outer loop by ktrlink and fmincon. This step is needed because ktrlink and fmincon require an objective function that only depends on the arguments of the minimization.
| DynamicDiscreteChoice.m | 147 | objfun = @(par)NegLogLikCCP(A,iX,calX,Pi,par(1:2),[delta(1);par(3)],rho,@flowpayoffs,@Psi,@FixedPointPsi,tFP) |
We are now ready to run the NFXP maximum likelihood procedure. If useKnitro equals 1, ktrlink will be called; otherwise, the Matlab procedure fmincon will be used.
| DynamicDiscreteChoice.m | 151 | if useKnitro |
| DynamicDiscreteChoice.m | 152 | options = optimset('Display','iter','GradObj','off','AlwaysHonorConstraints', ... |
| DynamicDiscreteChoice.m | 153 | 'bounds','TolFun',tc,'TolX',tp); |
| DynamicDiscreteChoice.m | 154 | [parMLE,fval,flag] = ktrlink(objfun,startvalues,[],[],[],[],lb,[],[],options); |
| DynamicDiscreteChoice.m | 155 | else |
| DynamicDiscreteChoice.m | 156 | options = optimset('Algorithm','interior-point','Display','iter','GradObj', ... |
| DynamicDiscreteChoice.m | 157 | 'on','AlwaysHonorConstraints','bounds','TolFun',tc,'TolX',tp,'DerivativeCheck','on'); |
| DynamicDiscreteChoice.m | 158 | [parMLE,fval,flag] = fmincon(objfun,startvalues,[],[],[],[],lb,[],[],options); |
| DynamicDiscreteChoice.m | 159 | end |
This gives maximum partial likelihood estimates of
. To calculate standard errors, we call NegLogLikCCP once more to estimate the corresponding Fisher information matrix and store this in gradgrad. Its inverse is an estimate of the maximum likelihood estimator's asymptotic variance-covariance matrix.
| DynamicDiscreteChoice.m | 164 | [nll grad gradgrad] = NegLogLikCCP(A,iX,calX,Pi,parMLE(1:2),[delta(1);parMLE(3)],rho,@flowpayoffs,@Psi,@FixedPointPsi,tFP); |
| DynamicDiscreteChoice.m | 165 | Vcov=inv(gradgrad) |
| DynamicDiscreteChoice.m | 166 | ses=sqrt(diag(Vcov)) |
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.m | 170 | disp([[beta;delta(2)] startvalues parMLE ses]); |
Extension to an Unknown Markov Transition Matrix for the State
Finally, consider the more realistic case that 
is not known. In this case, Rust (1994) suggests a two-stage procedure. In the first stage, we estimate 
using estimatePi and store the results in a 
matrix Pihat.
| DynamicDiscreteChoice.m | 175 | Pihat=estimatePi(iX,K) |
In the second stage, maximum partial likelihood estimates of 
can be computed using the NFXP procedure, with Pihat replacing Pi. This, with the question how the first-stage sampling error affects the precision of the second-stage estimator of 
, is left for the exercises.
6FStudent Exercises
Excercises to Work on
The 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 work with numerical derivatives by setting the option "GradObj" in fmincon or ktrlink to "off"; this would allow you to find the estimates, but not to compute their standard errors. Once you have coded up the analytical derivatives, you are advised to verify them against numerical derivatives, for example by setting the option "DerivativeCheck" in fmincon or ktrlink to "on".
- 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. Keep track of the estimates and report statistics such as their means and standard deviations. Note that you can also use this Monte Carlo setup in later questions to study the finite-sample performance of the various procedures studied.
- Implement the MPEC approach to maximum likelihood estimation of our structural model, as advocated by Judd and Su (2008).
- First, extend NegLogLikCCP so that it can be used to compute the log (partial) likelihood as a function of not only the parameters of interest but also the (finite number of) values of
. For example, allow for optional arguments capU0 and capU1 and, if these arguments are not empty, directly compute the choice probabilities from the implied values of 
, rather than computing 
by solving the model. - 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
, subject to the
constraint that 
(which can be specified using the
function Psi).
- First, extend NegLogLikCCP so that it can be used to compute the log (partial) likelihood as a function of not only the parameters of interest but also the (finite number of) values of
- 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
over the support of 
(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
by inverting the estimated choice probabilities (as proposed by Hotz and Miller (1993) and Hotz et al. (1994)).
- Extend NegLogLikCCP (or add functions) so that it optionally takes nonparametric estimates of
as inputs, computes the corresponding estimates of the entry and exit rules, and uses these estimates and the model to forwardly simulate 
and 
, and then 
, for each point 
in the sample. As the objective to be minimized, either specify a weighted distance between the nonparametric estimates of 
and the simulated values of 
, or minus a log pseudo-likelihood based on the choice probabilities implied by the simulated values of 
.
- 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).
- Add a new function for nonparametrically estimating the choice probabilities
- Extend the code for simulation and 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
and the other with entry cost 
. The firms are distributed across both unobserved entry cost types independently from all other variables in the model. We would like to estimate 
, and the share of agents with entry cost 
. To this end:
- Extend SimulateDP so that it randomly draws an entry cost from a distribution with two points of support,
and 
, and simulates choice data corresponding to the entry cost drawn, for each firm, independently across firms.
- Extend NegLokLikCCP 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.
- Extend SimulateDP so that it randomly draws an entry cost from a distribution with two points of support,
Additional Exercises
- Prove that
is a contraction.
- For the case of estimating the demand for differentiated products using NFXP maximum likelihood, as in Berry, Levinsohn, and Pakes (1995), Dub\'e, Fox, and Su (2008) 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 tFP on the one hand and tc and tp on the other hand to investigate this issue in the context of this package's firm entry and exit problem.
- Show the equivalence of policy iteration and the Newton-Kantorovich method for solving
for 
(see Appendix 1A.3). Extend FixedPointPsi so that it has the option to compute 
using the Newton-Kantorovich method (or a hybrid method, 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 
computed in an intermediate step of NegLogLikCCP.
- Compare the asymptotic standard errors to Monte Carlo standard errors.
- Compute the two-stage estimates of
for the case that 
is not known. Are the estimators of the standard errors that assume 
known consistent if 
is estimated instead? If not, extend the code so that it computes consistent standard error estimates for the case that 
is estimated. See Rust (1994).
- Is
identified in the example model? What about 
? To the extent that these are identified, extend your numerical experiments to include their estimation. See Abbring (2010) for results and references.
- Extend flowpayoffs.m and Psi.m so that they allow for models with different functional forms and a larger (but finite) number of choices and experiment with estimating these models.
1AContractions 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.
1A.1Definitions
A metric space is a set 
with a metric 
such that, for any 
,
-
if and only if 
,
-
, and
-
(triangle inequality).
A Cauchy sequence in 
is a sequence 
in 
such that, for each 
, there exists some 
such that 
for all 
.
A metric space 
is complete if every Cauchy sequence in 
converges to a point in 
.
A map 
is a contraction with modulus 
if 
for all 
.
1A.2Contraction Mapping (Banach Fixed Point) Theorem
If 
is a complete metric space and 
is a contraction, then there exists a unique 
such that 
. Furthermore, for any 
, the sequence 
with 
, 
, converges to 
.
Sketch of Proof
Suppose 
are such that 
and 
. Then, 
for some 
. Consequently, 
and 
.
Because 
is a contraction, any 
it generates is a Cauchy sequence. Because 
is complete, it converges to some 
. Because contractions are continuous, 
.
1A.3Computing the Fixed Point of a Contraction
The method of successive approximations directly applies the Contraction Mapping Theorem and
approximates the fixed point 
of a contraction 
with the sequence 
generated from 
on a finitely discretized
state space. This method has global but linear convergence.
The Newton-Kantorovich method searches a zero of 
; it approximates 
with the sequence generated from
(0)
with 
the Fr\'echet derivative of 
at 
(with a discrete state space, this is simply the
linear map defined by the appropriate finite dimensional Jacobian, as in NegLogLikCCP).
1A.4Blackwell's Sufficient Conditions for a Contraction
Let 
be the space of functions 
such that 
( 
bounded), with metric 
. Suppose that, for some 
, 
satisfies
- (monotonicity)
and
- (discounting)
for all 
such that 
and all 
. Then, 
is a contraction with modulus 
.
Sketch of Proof
Let 
be such that 
. Note that 
, so that 
by monot. & discounting. Similarly, 
. Thus, 
.
1A.5Application to Dynamic Programming
The choice-specific value function 
is a fixed point of 
, which satisfies Blackwell's sufficient conditions with 
equal to the discount factor.
- Under regularity conditions such that
is bounded, this ensures that it is the unique fixed point of a contraction on the complete (Banach) space 
of bounded functions.
- If
is closed and 
for all 
, then 
is a contraction on the complete subspace 
, so that 
(examples: monotonicity, continuity, etcetera).
-
can be computed by value iteration (successive approximations), policy iteration (the Newton-Kantorovich Method), or a hybrid algorithm as in Rust (1987).
2BMiscellaneous Utilities
This section documents some miscellaneous utilities that are called by the main functions. At this point, it only includes RandomDiscr.
The function RandomDiscr returns a random draw from the distribution of 
; with 
independently distributed with (not necessarily identical) disctributions on 
for some 
.
| RandomDiscr.m | 7 | function iX = RandomDiscr(P) |
It requires one input argument:
- P
- a
matrix with 
th element 
.
It returns
- iX
- a
vector with a random draw 
from 
.
The function RandomDiscr first stores the number 
of support points in a scalar K and the number 
of random variables in a scalar N.
| RandomDiscr.m | 19 | K = size(P,1); |
| RandomDiscr.m | 20 | N = size(P,2); |
Then, it creates a 
matrix ru with 
identical rows, containing 
independent draws from a standard uniform distribution, and computes the 
matrix cumP with 
th element 
.
| RandomDiscr.m | 24 | ru = ones(K-1,1)*random('unif',zeros(1,N),ones(1,N)); |
| RandomDiscr.m | 25 | cumP = cumsum(P); |
Finally, for each 
, it sets 
equal to 1 plus the number of elements of 
} that are weakly smaller than the uniform random draw in the 
th column of ru.
| RandomDiscr.m | 29 | iX = sum([ones(1,N);cumP(1:K-1,:)<=ru]); |
- Abbring, Jaap H. (2010). Identification of dynamic discrete choice models. Annual Review of Economics 2, 367-394.
- 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.
- Campbell, Jeffrey R. and Ryan Peters (2010). Comments ++.
- Dub\'e, Jean-Pierre, Jeremy T. Fox, and Che-Lin Su (2008). Improving the numerical performance of BLP structural demand estimators. mimeo. University of Chicago.
- Dixit, Avinash K. (1989). Entry and exit decisions under uncertainty. Journal of Political Economy 97(3), 620-638.
- 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.
- Judd, Kenneth L. 1998. Numerical Methods in Economics. MIT Press. Cambridge, MA.
- Judd, Kenneth L. and Che-Lin Su (2008). Constrained optimization approaches to estimation of structural models. mimeo. University of Chicago.
- Hotz, Joseph V. and David A. Miller (1993). Conditional choice probabilities and the estimation of dynamic models. Review of Economic Studies 60, 497-529.
- Hotz, Joseph V., David A. Miller, S. Sanders and J. Smith (1994). A simulation estimator for dynamic models of discrete choice. Review of Economic Studies 61(2), 265-289.
- Keane, Michael P. and Kenneth I. Wolpin (1997). The career decisions of young men. Journal of Political Economy 105(3), 473-522.
- McFadden, Daniel (1981). Econometric models of probabilistic choice. In C. Manski and D. McFadden (Eds.). Discrete Data with Econometric Applications. MIT Press. Cambridge, MA.
- 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.
