Dynamic Structural Econometrics - Unobserved heterogeneity and CCP estimation

The Expectation-Maximization (EM) Algorithm

rng(2023);

nobs=1000;

p=rand([nobs,1]);

Y=(p>0.6).*poissrnd(2,[nobs,1]) + (p<=0.6).*poissrnd(8,[nobs,1]);

histogram(Y,'Normalization',"probability")

poiss_hat = poissfit(Y);

nbin_hat = nbinfit(Y);

x=0:16;

histogram(Y,'Normalization','pdf','FaceColor',[0.4,0.4,0.4]);

hold on

plot(x,poisspdf(x,poiss_hat));

plot(x,nbinpdf(x,nbin_hat(1),nbin_hat(2)),'--');

hold off;

legend('Empirical distribution','Poisson estimates','Neg. Bin. estimates','location','northeast')

A common way to estimate these distributions is to consider a mixture model. That is, we will assume that each draw of Y is obtained from distribution with probability and distribution with probability . The mixing parameters govern the distribution over types, while the distributions describe the distribution of outcomes within each type. Typically, we will parameterize the distributions and with parameters . Noting the restriction that , we have three parameters we need to estimate: , , and .

The log-likelihood for observation can be expressed as

Practical applications have found that direct maximum likelihood approaches to estimating these parameters perform poorly -- convergence may be slow, there can be many local optima, and the estimates may be extremely sensitive to starting values for the parameters. The alternative approach of Dempster et al (1977) is to update the parameters π and in separate stages. The main features of the estimator can be described in four steps:

1. Fix initial values for , and and set convergence criteria.
2. Given the current values for , compute the probability of each observation type, i.e., calculate Prfor . This yields an updated value for π and pointwise estimates for the probability that each observation i is of type j.
3. Hold fixed the type probabilities for each computed in step 2 and maximize the likelihood with respect to given these probabilities. This yields updated values for .
4. Repeat steps 2-3 until convergence (per-iteration change in the overall likelihood is "small").

We will implement this algorithm on the data set above using the assumptions that (i) there are two components in the mixture (i.e., observations are each associated with one of two unobserved types) and (ii) the two types each follow a Poisson distribution with parameter .

Expectation step.

The expectation step takes the current values of as fixed and calculates the probability of observation being of type 1 () or type 2 (). We can write the density function of Y conditional on θ, π, and the (unobserved) as

Similarly, the distribution of conditional on θ and π is

where . These equations together with Bayes theorem imply that the conditional distribution of the are the proportional likelihood of each mixture component j weighted by the current estimates of the mixture weights λ:

In the code above, the "temp_pdf" variable creates an matrix of likelihood values, where the first column is the set of likelihoods associated with and the second column is the set of likelihoods associated with :

Maximization step.

The maximization step performs a standard maximum likelihood procedure in which the probabilities of being in each group are treated as known and fixed.

Because the likelihood is additively separable in θ and λ, the maximization step can be done in two parts. Updating π is simple: maximizing with respect to π simply requires taking the average of the probability of being each type:

That is, we simply take the mean of each column of the "EZ" matrix we computed above.

The last step in each iteration is to calculate the overall likelihood and check for convergence. If is less than the convergence criteria, we stop. Otherwise, we proceed to the next iteration.

x=0:20;

Y_fit=PiHat(1)*poisspdf(x,ThetaHat(1))+PiHat(2)*poisspdf(x,ThetaHat(2));

plot(x,Y_fit)

hold on

histogram(Y,'Normalization','probability','FaceColor',[0.4,0.4,0.4]);

hold off

legend('Estimated mixture distribution','Empirical distribution','location','northeast')

Using the EM algorithm in a dynamic discrete choice model

Model description and setup

• Time is discrete, with denoting each period, with .
• Each period, Harold Zurcher chooses between two actions. He can either replace the bus engine, , or not replace . These are mutually exclusive, so that .
• Zurcher's payoff depends on a single state variable that represents the mileage of the bus. Mileage advances by a random amount if the bus engine is not replaced and is set to zero if the bus engine is replaced. Denote by the distribution over next period mileage conditional on current mileage after making choice d. When the bus engine is replaced, mileage resets to zero, so that for all . Mileage accumulation when the engine is not replaced follows a discrete analog of the exponential distribution up to some maximum mileage where mileage increment is .

• In addition to mileage , there is a permanent and unobservable bus-specific type . Assume that the deterministic component of flow payoffs for each action for a bus of type s are given by

• Payoffs also have an additively separable idiosyncratic component so that the total flow payoff from choice j in period t is . We assume that are i.i.d draws from a Type I EV distribution.
• Denoting the optimal decision rule at t by , we define the ex-ante value function for state in period t as

• The conditional value function is the value function for choice j in period t and state less the utility shock . In the case of Zurcher's bus engine replacement problem, these functions can be expressed as

• Note that all t subscripts have been removed from value functions and state transitions: this is a stationary, infinite horizon problem.
• In the binary choice setting with Type I EV preference shocks we can express conditional choice probabilities very simply:

• Representation: finally, we can use the results of Hotz and Miller (1993) to express the difference between the value function and conditional value function :

where when Type I EV, where is the Euler-Mascheroni constant.

Estimation

We begin by setting initial values for the parameters and , and use the initial values for to initialize the conditional choice probabilities .

Expectation step.

With initial values for the parameters established, we first need to calculate the likelihood for each observation bus n, time period t, and bus type s. Recall that in the finite Poisson mixture in the example above, we updated the guesses for the probability that observation i was type j by using the π-weighted likelihoods of observing outcome conditional on being type j and having type-specific parameters . With the dynamic bus engine replacement problem, we proceed similarly. However, we have to remember that the relevant likelihood for establishing whether bus n is of type s is the sequence of replacement decisions made for that bus. Therefore, we calculate the individual likelihood components for each period, bus, and type, and then we take the product across time periods t to compute the overall likelihood associated with each bus type s.

The likelihood associated with each bus n conditional on engine type s is the product of the likelihood of the sequence of decisions conditional on the states , engine type s, conditional choice probabilities , and utility parameters θ: . This is the term that captures the observation-level likelihood terms from our simple EM algorithm example above. As before, we weight the engine-type specific likelihoods by the mixing parameter yielding an overal per-bus estimate of the probability of being engine type s of

Maximization step.

To update the utility parameters θ, we first need to update the conditional choice probabilities since we are using CCP methods to calculate the likelihood terms. We update the conditional choice probabilities using the bin estimator after applying weights via the per-bus probability of being in state s:

The last maximization step is to update the utility parameters θ. Here we maximize the likelihood, again weighting each likelihood component by the probability that bus n is of type s using the estimates for obtained in the previous step:

Important disclaimer: in the data generation function, we have assigned exogenous state variables in the initial period . This is unlikely to hold in general. Suppose all buses had started with zero mileage in period 0 but we, the econometricians, only observed the buses starting in period 10. Because the mileage accumulation and replacement decisions will be correlated with the unobserved type s, we will have an initial conditions problem: the first state that we observed the bus, , will be correlated with the unobserved type, and we will need to control for this in order to obtain consistent estimates. There are several ways to account for this; see the supplementary material for Arcidiacono and Miller (2011) for one such implementation.

Functions

References