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

Jaap H. Abbring§Tobias J. Klein†

This version: June 7, 2011

###

This document supports the first Matlab computing sessions in our PhD elective course Empirical Industrial Organization II in CentER Tilburg's Research Master in Economics program (230323) and in the Finnish Doctoral Programme in Economics (SEIO11). It contains some notes on the theory of dynamic discrete choice models and on methods for their computation and estimation. It is centered around some basic Matlab code for solving, simulating, and empirically analyzing a simple dynamic discrete choice model. Student exercises ask students to extend this code to apply different and more advanced computational and econometric methods to a wider range of models.

✶ | Thanks to Jeffrey R. Campbell for text on and help with Comments++ and to Nan Yang for comments.
§ | CentER, Department of Econometrics & OR, Tilburg University, P.O. Box 90153, 5000 LE Tilburg, The Netherlands. E-mail: j.h.Abbring@uvt.nl. http://jabbring.home.xs4all.nl.
† | CentER, Department of Econometrics & OR, Tilburg University.

###

- 1 A A Firm Entry and Exit Problem
  - 1 A.1 Flow Payoffs
  - 1 A.2 Bellman Operator
- 2 B Solving the Decision Problem
- 3 C Data Simulation
- 4 D Estimation
  - 4 D.1 First Stage: State Transitions
  - 4 D.2 Second Stage: Choices
- 5 E The Script that Puts It All Together
- 6 F Student Exercises

- 1 A Contractions and Dynamic Programming
  - 1 A.1 Definitions
  - 1 A.2 Contraction Mapping (Banach Fixed Point) Theorem
  - 1 A.3 Computing the Fixed Point of a Contraction
  - 1 A.4 Blackwell's Sufficient Conditions for a Contraction
  - 1 A.5 Application to Dynamic Programming
- 2 B Miscellaneous Utilities

#### 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 2 B 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 3 C documents how SimulateDP can be used to simulate data for the decision problem set up by flowpayoffs and Psi. Section 4 D documents functions for estimating the transition matrix and for computing the log partial likelihood for conditional choice data. Section 5 E 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 F contains a range of student exercises. Appendix 1 A provides mathematical background by discussing some theory of contractions. Appendix 2 B documents procedures that are called by the main routines, but that are of limited substantial interest, such as RandomDiscr.

# 1 A A 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.

## 1 A.1 Flow 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.

## 1 A.2 Bellman 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).

# 2 B Solving 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 1 A.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.

# 3 C Data 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 2 B) 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

# 4 D Estimation

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.

## 4 D.1 First 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.

## 4 D.2 Second 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 5 E'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 .

# 5 E The 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 3 C 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.

# 6 F Student 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".

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. 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.
2. 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).
   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?
3. 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).
   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)).
4. 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.

#### Additional Exercises

1. Prove that is a contraction.
2. 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.
3. Show the equivalence of policy iteration and the Newton-Kantorovich method for solving for (see Appendix 1 A.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.
4. Compare the asymptotic standard errors to Monte Carlo standard errors.
5. 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).
6. 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.
7. 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.

###

# 1 A Contractions 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.

## 1 A.1 Definitions

A metric space is a set with a metric such that, for any ,

1. if and only if ,
2. , and
3. (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 .

## 1 A.2 Contraction 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, .

## 1 A.3 Computing 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).

## 1 A.4 Blackwell's Sufficient Conditions for a Contraction

Let be the space of functions such that ( bounded), with metric . Suppose that, for some , satisfies

1. (monotonicity) and
2. (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, .

## 1 A.5 Application 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).

# 2 B Miscellaneous 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.

###

1 | A similar option is available on Microsoft Internet Explorer. We use Cascading Style Sheets (CSS2) for all formatting. Since IE implements this standard poorly, we cannot guarantee that it will display this document well.
2 | We have successfully printed this using Firefox for the Macintosh and Firefox for Linux. Sadly, not all browsers are so cooperative. This document apparently triggers a bug in Safari that causes multiple blank pages to be printed after each actual page.
3 | Please see the manual for the Matlab-Knitro function ktrlink or the Matlab function fmincon to learn about the tolerances and other input arguments to the minimization function used in the outer loop of the NFXP procedure.
4 | Note that both ktrlink and fmincon 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 for to ensure that ). The functions ktrlink and 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 (2008) (see Section 6 F).

DynamicDiscreteChoice.m | 93 | T = 100
DynamicDiscreteChoice.m | 94 | N = 1000
DynamicDiscreteChoice.m | 98 | tFP = 1e-10
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
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;
DynamicDiscreteChoice.m | 126 | [A,iX] = SimulateDP(DeltaU,Pi,T,N);
DynamicDiscreteChoice.m | 132 | useKnitro = 0;
DynamicDiscreteChoice.m | 133 | tc = 1E-6
DynamicDiscreteChoice.m | 134 | tp = 1E-6
DynamicDiscreteChoice.m | 138 | startvalues = [-1;-0.1;0.5];
DynamicDiscreteChoice.m | 142 | lb = -Inf*ones(size(startvalues));
DynamicDiscreteChoice.m | 143 | lb(3) = 0;
DynamicDiscreteChoice.m | 147 | objfun = @(par)NegLogLikCCP(A,iX,calX,Pi,par(1:2),[delta(1);par(3)],rho,@flowpayoffs,@Psi,@FixedPointPsi,tFP)
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
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))
DynamicDiscreteChoice.m | 170 | disp([[beta;delta(2)] startvalues parMLE ses]);
DynamicDiscreteChoice.m | 175 | Pihat=estimatePi(iX,K)
flowpayoffs.m | 7 | function [u0,u1] = flowpayoffs(calX,beta,delta)
flowpayoffs.m | 24 | K = size(calX,1);
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];
Psi.m | 15 | function [capU0,capU1] = Psi(capU0,capU1,u0,u1,Pi,rho)
Psi.m | 33 | R0 = log(exp(capU0(:,1))+exp(capU1(:,1)));
Psi.m | 34 | R1 = log(exp(capU0(:,2))+exp(capU1(:,2)));
Psi.m | 38 | capU0 = u0 + rho*Pi*R0*[1 1];
Psi.m | 39 | capU1 = u1 + rho*Pi*R1*[1 1];
FixedPointPsi.m | 8 | function [capU0,capU1] = FixedPointPsi(u0,u1,Pi,rho,tFP,Psi,capU0,capU1)
FixedPointPsi.m | 28 | K = size(Pi,1);
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
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
SimulateDP.m | 7 | function [A,iX] = SimulateDP(DeltaU,Pi,T,N)
SimulateDP.m | 23 | K = size(Pi,1);
SimulateDP.m | 39 | oneminPi = eye(K)-Pi';
SimulateDP.m | 40 | Pinf = inv([oneminPi(1:K-1,:);ones(1,K)])*[zeros(K-1,1);1];
SimulateDP.m | 44 | iX = RandomDiscr(Pinf*ones(1,N));
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;
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
estimatePi.m | 7 | function Pihat = estimatePi(iX,K)
estimatePi.m | 20 | T=size(iX,1);
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
NegLogLikCCP.m | 7 | function [nll grad gradgrad] = NegLogLikCCP(A,iX,calX,Pi,beta,delta,rho,flowpayoffs,Psi,FixedPointPsi,tFP)
NegLogLikCCP.m | 34 | K = size(calX,1);
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));
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)));
NegLogLikCCP.m | 56 | if nargout>=2
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]];
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)]];
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,:);
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
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
RandomDiscr.m | 7 | function iX = RandomDiscr(P)
RandomDiscr.m | 19 | K = size(P,1);
RandomDiscr.m | 20 | N = size(P,2);
RandomDiscr.m | 24 | ru = ones(K-1,1)*random('unif',zeros(1,N),ones(1,N));
RandomDiscr.m | 25 | cumP = cumsum(P);
RandomDiscr.m | 29 | iX = sum([ones(1,N);cumP(1:K-1,:)<=ru]);