Maximum Likelihood Estimation¶

Contents¶

• Maximum Likelihood Estimation
  • Overview
  • Set Up and Assumptions
  • Conditional Distributions
  • Maximum Likelihood Estimation
  • MLE with Numerical Methods
  • Maximum Likelihood Estimation with statsmodels
  • Summary
  • Exercises
  • Solutions

Overview¶

In a previous lecture, we estimated the relationship between dependent and explanatory variables using linear regression.

But what if a linear relationship is not an appropriate assumption for our model?

One widely used alternative is maximum likelihood estimation, which involves specifying a class of distributions, indexed by unknown parameters, and then using the data to pin down these parameter values.

The benefit relative to linear regression is that it allows more flexibility in the probabilistic relationships between variables.

Here we illustrate maximum likelihood by replicating Daniel Treisman’s (2016) paper, Russia’s Billionaires, which connects the number of billionaires in a country to its economic characteristics.

The paper concludes that Russia has a higher number of billionaires than economic factors such as market size and tax rate predict.

We’ll require the following imports:

Prerequisites¶

We assume familiarity with basic probability and multivariate calculus.

Set Up and Assumptions¶

Let’s consider the steps we need to go through in maximum likelihood estimation and how they pertain to this study.

Flow of Ideas¶

The first step with maximum likelihood estimation is to choose the probability distribution believed to be generating the data.

More precisely, we need to make an assumption as to which parametric class of distributions is generating the data.

• e.g., the class of all normal distributions, or the class of all gamma distributions.

Each such class is a family of distributions indexed by a finite number of parameters.

• e.g., the class of normal distributions is a family of distributions indexed by its mean $ \mu \in (-\infty, \infty) $ and standard deviation $ \sigma \in (0, \infty) $.

We’ll let the data pick out a particular element of the class by pinning down the parameters.

The parameter estimates so produced will be called maximum likelihood estimates.

Counting Billionaires¶

Treisman [Tre16] is interested in estimating the number of billionaires in different countries.

The number of billionaires is integer-valued.

Hence we consider distributions that take values only in the nonnegative integers.

(This is one reason least squares regression is not the best tool for the present problem, since the dependent variable in linear regression is not restricted to integer values)

One integer distribution is the Poisson distribution, the probability mass function (pmf) of which is

$$ f(y) = \frac{\mu^{y}}{y!} e^{-\mu}, \qquad y = 0, 1, 2, \ldots, \infty $$

We can plot the Poisson distribution over $ y $ for different values of $ \mu $ as follows

Notice that the Poisson distribution begins to resemble a normal distribution as the mean of $ y $ increases.

Let’s have a look at the distribution of the data we’ll be working with in this lecture.

Treisman’s main source of data is Forbes’ annual rankings of billionaires and their estimated net worth.

The dataset mle/fp.dta can be downloaded here or from its AER page.

Using a histogram, we can view the distribution of the number of billionaires per country, numbil0, in 2008 (the United States is dropped for plotting purposes)

From the histogram, it appears that the Poisson assumption is not unreasonable (albeit with a very low $ \mu $ and some outliers).

Conditional Distributions¶

In Treisman’s paper, the dependent variable — the number of billionaires $ y_i $ in country $ i $ — is modeled as a function of GDP per capita, population size, and years membership in GATT and WTO.

Hence, the distribution of $ y_i $ needs to be conditioned on the vector of explanatory variables $ \mathbf{x}_i $.

The standard formulation — the so-called poisson regression model — is as follows:

$$ f(y_i \mid \mathbf{x}_i) = \frac{\mu_i^{y_i}}{y_i!} e^{-\mu_i}; \qquad y_i = 0, 1, 2, \ldots , \infty . \tag{1} $$

$$ \text{where}\ \mu_i = \exp(\mathbf{x}_i' \boldsymbol{\beta}) = \exp(\beta_0 + \beta_1 x_{i1} + \ldots + \beta_k x_{ik}) $$

To illustrate the idea that the distribution of $ y_i $ depends on $ \mathbf{x}_i $ let’s run a simple simulation.

We use our poisson_pmf function from above and arbitrary values for $ \boldsymbol{\beta} $ and $ \mathbf{x}_i $

We can see that the distribution of $ y_i $ is conditional on $ \mathbf{x}_i $ ($ \mu_i $ is no longer constant).

Maximum Likelihood Estimation¶

In our model for number of billionaires, the conditional distribution contains 4 ($ k = 4 $) parameters that we need to estimate.

We will label our entire parameter vector as $ \boldsymbol{\beta} $ where

$$ \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{bmatrix} $$

To estimate the model using MLE, we want to maximize the likelihood that our estimate $ \hat{\boldsymbol{\beta}} $ is the true parameter $ \boldsymbol{\beta} $.

Intuitively, we want to find the $ \hat{\boldsymbol{\beta}} $ that best fits our data.

First, we need to construct the likelihood function $ \mathcal{L}(\boldsymbol{\beta}) $, which is similar to a joint probability density function.

Assume we have some data $ y_i = \{y_1, y_2\} $ and $ y_i \sim f(y_i) $.

If $ y_1 $ and $ y_2 $ are independent, the joint pmf of these data is $ f(y_1, y_2) = f(y_1) \cdot f(y_2) $.

If $ y_i $ follows a Poisson distribution with $ \lambda = 7 $, we can visualize the joint pmf like so

Similarly, the joint pmf of our data (which is distributed as a conditional Poisson distribution) can be written as

$$ f(y_1, y_2, \ldots, y_n \mid \mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n; \boldsymbol{\beta}) = \prod_{i=1}^{n} \frac{\mu_i^{y_i}}{y_i!} e^{-\mu_i} $$

$ y_i $ is conditional on both the values of $ \mathbf{x}_i $ and the parameters $ \boldsymbol{\beta} $.

The likelihood function is the same as the joint pmf, but treats the parameter $ \boldsymbol{\beta} $ as a random variable and takes the observations $ (y_i, \mathbf{x}_i) $ as given

$$ \begin{split} \mathcal{L}(\beta \mid y_1, y_2, \ldots, y_n \ ; \ \mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n) = & \prod_{i=1}^{n} \frac{\mu_i^{y_i}}{y_i!} e^{-\mu_i} \\ = & f(y_1, y_2, \ldots, y_n \mid \ \mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n ; \beta) \end{split} $$

Now that we have our likelihood function, we want to find the $ \hat{\boldsymbol{\beta}} $ that yields the maximum likelihood value

$$ \underset{\boldsymbol{\beta}}{\max} \mathcal{L}(\boldsymbol{\beta}) $$

In doing so it is generally easier to maximize the log-likelihood (consider differentiating $ f(x) = x \exp(x) $ vs. $ f(x) = \log(x) + x $).

Given that taking a logarithm is a monotone increasing transformation, a maximizer of the likelihood function will also be a maximizer of the log-likelihood function.

In our case the log-likelihood is

$$ \begin{split} \log{ \mathcal{L}} (\boldsymbol{\beta}) = \ & \log \Big( f(y_1 ; \boldsymbol{\beta}) \cdot f(y_2 ; \boldsymbol{\beta}) \cdot \ldots \cdot f(y_n ; \boldsymbol{\beta}) \Big) \\ = & \sum_{i=1}^{n} \log{f(y_i ; \boldsymbol{\beta})} \\ = & \sum_{i=1}^{n} \log \Big( {\frac{\mu_i^{y_i}}{y_i!} e^{-\mu_i}} \Big) \\ = & \sum_{i=1}^{n} y_i \log{\mu_i} - \sum_{i=1}^{n} \mu_i - \sum_{i=1}^{n} \log y! \end{split} $$

The MLE of the Poisson to the Poisson for $ \hat{\beta} $ can be obtained by solving

$$ \underset{\beta}{\max} \Big( \sum_{i=1}^{n} y_i \log{\mu_i} - \sum_{i=1}^{n} \mu_i - \sum_{i=1}^{n} \log y! \Big) $$

However, no analytical solution exists to the above problem – to find the MLE we need to use numerical methods.

MLE with Numerical Methods¶

Many distributions do not have nice, analytical solutions and therefore require numerical methods to solve for parameter estimates.

One such numerical method is the Newton-Raphson algorithm.

Our goal is to find the maximum likelihood estimate $ \hat{\boldsymbol{\beta}} $.

At $ \hat{\boldsymbol{\beta}} $, the first derivative of the log-likelihood function will be equal to 0.

Let’s illustrate this by supposing

$$ \log \mathcal{L(\beta)} = - (\beta - 10) ^2 - 10 $$

The plot shows that the maximum likelihood value (the top plot) occurs when $ \frac{d \log \mathcal{L(\boldsymbol{\beta})}}{d \boldsymbol{\beta}} = 0 $ (the bottom plot).

Therefore, the likelihood is maximized when $ \beta = 10 $.

We can also ensure that this value is a maximum (as opposed to a minimum) by checking that the second derivative (slope of the bottom plot) is negative.

The Newton-Raphson algorithm finds a point where the first derivative is 0.

To use the algorithm, we take an initial guess at the maximum value, $ \beta_0 $ (the OLS parameter estimates might be a reasonable guess), then

1. Use the updating rule to iterate the algorithm
   $$ \boldsymbol{\beta}_{(k+1)} = \boldsymbol{\beta}_{(k)} - H^{-1}(\boldsymbol{\beta}_{(k)})G(\boldsymbol{\beta}_{(k)}) $$ where:
   $$ \begin{aligned} G(\boldsymbol{\beta}_{(k)}) = \frac{d \log \mathcal{L(\boldsymbol{\beta}_{(k)})}}{d \boldsymbol{\beta}_{(k)}} \\ H(\boldsymbol{\beta}_{(k)}) = \frac{d^2 \log \mathcal{L(\boldsymbol{\beta}_{(k)})}}{d \boldsymbol{\beta}_{(k)}d \boldsymbol{\beta}'_{(k)}} \end{aligned} $$
2. Check whether $ \boldsymbol{\beta}_{(k+1)} - \boldsymbol{\beta}_{(k)} < tol $
   • If true, then stop iterating and set $ \hat{\boldsymbol{\beta}} = \boldsymbol{\beta}_{(k+1)} $
   • If false, then update $ \boldsymbol{\beta}_{(k+1)} $

As can be seen from the updating equation, $ \boldsymbol{\beta}_{(k+1)} = \boldsymbol{\beta}_{(k)} $ only when $ G(\boldsymbol{\beta}_{(k)}) = 0 $ ie. where the first derivative is equal to 0.

(In practice, we stop iterating when the difference is below a small tolerance threshold)

Let’s have a go at implementing the Newton-Raphson algorithm.

First, we’ll create a class called PoissonRegression so we can easily recompute the values of the log likelihood, gradient and Hessian for every iteration

Our function newton_raphson will take a PoissonRegression object that has an initial guess of the parameter vector $ \boldsymbol{\beta}_0 $.

The algorithm will update the parameter vector according to the updating rule, and recalculate the gradient and Hessian matrices at the new parameter estimates.

Iteration will end when either:

• The difference between the parameter and the updated parameter is below a tolerance level.
• The maximum number of iterations has been achieved (meaning convergence is not achieved).

So we can get an idea of what’s going on while the algorithm is running, an option display=True is added to print out values at each iteration.

Let’s try out our algorithm with a small dataset of 5 observations and 3 variables in $ \mathbf{X} $.

As this was a simple model with few observations, the algorithm achieved convergence in only 6 iterations.

You can see that with each iteration, the log-likelihood value increased.

Remember, our objective was to maximize the log-likelihood function, which the algorithm has worked to achieve.

Also, note that the increase in $ \log \mathcal{L}(\boldsymbol{\beta}_{(k)}) $ becomes smaller with each iteration.

This is because the gradient is approaching 0 as we reach the maximum, and therefore the numerator in our updating equation is becoming smaller.

The gradient vector should be close to 0 at $ \hat{\boldsymbol{\beta}} $

The iterative process can be visualized in the following diagram, where the maximum is found at $ \beta = 10 $

Note that our implementation of the Newton-Raphson algorithm is rather basic — for more robust implementations see, for example, scipy.optimize.

Maximum Likelihood Estimation with statsmodels¶

Now that we know what’s going on under the hood, we can apply MLE to an interesting application.

We’ll use the Poisson regression model in statsmodels to obtain a richer output with standard errors, test values, and more.

statsmodels uses the same algorithm as above to find the maximum likelihood estimates.

Before we begin, let’s re-estimate our simple model with statsmodels to confirm we obtain the same coefficients and log-likelihood value.

Now let’s replicate results from Daniel Treisman’s paper, Russia’s Billionaires, mentioned earlier in the lecture.

Treisman starts by estimating equation (1), where:

• $ y_i $ is $ {number\ of\ billionaires}_i $
• $ x_{i1} $ is $ \log{GDP\ per\ capita}_i $
• $ x_{i2} $ is $ \log{population}_i $
• $ x_{i3} $ is $ {years\ in\ GATT}_i $ – years membership in GATT and WTO (to proxy access to international markets)

The paper only considers the year 2008 for estimation.

We will set up our variables for estimation like so (you should have the data assigned to df from earlier in the lecture)

Then we can use the Poisson function from statsmodels to fit the model.

We’ll use robust standard errors as in the author’s paper

Success! The algorithm was able to achieve convergence in 9 iterations.

Our output indicates that GDP per capita, population, and years of membership in the General Agreement on Tariffs and Trade (GATT) are positively related to the number of billionaires a country has, as expected.

Let’s also estimate the author’s more full-featured models and display them in a single table

The output suggests that the frequency of billionaires is positively correlated with GDP per capita, population size, stock market capitalization, and negatively correlated with top marginal income tax rate.

To analyze our results by country, we can plot the difference between the predicted an actual values, then sort from highest to lowest and plot the first 15

As we can see, Russia has by far the highest number of billionaires in excess of what is predicted by the model (around 50 more than expected).

Treisman uses this empirical result to discuss possible reasons for Russia’s excess of billionaires, including the origination of wealth in Russia, the political climate, and the history of privatization in the years after the USSR.

Summary¶

In this lecture, we used Maximum Likelihood Estimation to estimate the parameters of a Poisson model.

statsmodels contains other built-in likelihood models such as Probit and Logit.

For further flexibility, statsmodels provides a way to specify the distribution manually using the GenericLikelihoodModel class - an example notebook can be found here.

Exercises¶

Exercise 1¶

Suppose we wanted to estimate the probability of an event $ y_i $ occurring, given some observations.

We could use a probit regression model, where the pmf of $ y_i $ is

$$ \begin{aligned} f(y_i; \boldsymbol{\beta}) = \mu_i^{y_i} (1-\mu_i)^{1-y_i}, \quad y_i = 0,1 \\ \text{where} \quad \mu_i = \Phi(\mathbf{x}_i' \boldsymbol{\beta}) \end{aligned} $$

$ \Phi $ represents the cumulative normal distribution and constrains the predicted $ y_i $ to be between 0 and 1 (as required for a probability).

$ \boldsymbol{\beta} $ is a vector of coefficients.

Following the example in the lecture, write a class to represent the Probit model.

To begin, find the log-likelihood function and derive the gradient and Hessian.

The scipy module stats.norm contains the functions needed to compute the cmf and pmf of the normal distribution.

Exercise 2¶

Use the following dataset and initial values of $ \boldsymbol{\beta} $ to estimate the MLE with the Newton-Raphson algorithm developed earlier in the lecture

$$ \mathbf{X} = \begin{bmatrix} 1 & 2 & 4 \\ 1 & 1 & 1 \\ 1 & 4 & 3 \\ 1 & 5 & 6 \\ 1 & 3 & 5 \end{bmatrix} \quad y = \begin{bmatrix} 1 \\ 0 \\ 1 \\ 1 \\ 0 \end{bmatrix} \quad \boldsymbol{\beta}_{(0)} = \begin{bmatrix} 0.1 \\ 0.1 \\ 0.1 \end{bmatrix} $$

Verify your results with statsmodels - you can import the Probit function with the following import statement

Note that the simple Newton-Raphson algorithm developed in this lecture is very sensitive to initial values, and therefore you may fail to achieve convergence with different starting values.

Solutions¶

Exercise 1¶

The log-likelihood can be written as

$$ \log \mathcal{L} = \sum_{i=1}^n \big[ y_i \log \Phi(\mathbf{x}_i' \boldsymbol{\beta}) + (1 - y_i) \log (1 - \Phi(\mathbf{x}_i' \boldsymbol{\beta})) \big] $$

Using the fundamental theorem of calculus, the derivative of a cumulative probability distribution is its marginal distribution

$$ \frac{ \partial} {\partial s} \Phi(s) = \phi(s) $$

where $ \phi $ is the marginal normal distribution.

The gradient vector of the Probit model is

$$ \frac {\partial \log \mathcal{L}} {\partial \boldsymbol{\beta}} = \sum_{i=1}^n \Big[ y_i \frac{\phi(\mathbf{x}'_i \boldsymbol{\beta})}{\Phi(\mathbf{x}'_i \boldsymbol{\beta)}} - (1 - y_i) \frac{\phi(\mathbf{x}'_i \boldsymbol{\beta)}}{1 - \Phi(\mathbf{x}'_i \boldsymbol{\beta)}} \Big] \mathbf{x}_i $$

The Hessian of the Probit model is

$$ \frac {\partial^2 \log \mathcal{L}} {\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}'} = - \sum_{i=1}^n \phi (\mathbf{x}_i' \boldsymbol{\beta}) \Big[ y_i \frac{ \phi (\mathbf{x}_i' \boldsymbol{\beta}) + \mathbf{x}_i' \boldsymbol{\beta} \Phi (\mathbf{x}_i' \boldsymbol{\beta}) } { [\Phi (\mathbf{x}_i' \boldsymbol{\beta})]^2 } + (1 - y_i) \frac{ \phi_i (\mathbf{x}_i' \boldsymbol{\beta}) - \mathbf{x}_i' \boldsymbol{\beta} (1 - \Phi (\mathbf{x}_i' \boldsymbol{\beta})) } { [1 - \Phi (\mathbf{x}_i' \boldsymbol{\beta})]^2 } \Big] \mathbf{x}_i \mathbf{x}_i' $$

Using these results, we can write a class for the Probit model as follows

Exercise 2¶