---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.10.3
kernelspec:
  display_name: Python 3
  language: python
  name: pymc
---

```{code-cell} ipython3
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import gamma
from scipy.stats import uniform
from scipy.optimize import minimize
import statsmodels.formula.api as smf 
import seaborn as sbn

np.random.seed(12578)
```

# Econometric routines and statistical distributions in Python

## Statsmodels

The statsmodels package has some of the functionality of stata.  It can estimate most models you would see as a student taking econometrics, time series, and cross section here at William and Mary.  Here is an example OLS regression:

```{code-cell} ipython3
# load data and create dataframe
tobias_koop=pd.read_csv('https://rlhick.people.wm.edu/econ407/data/tobias_koop_t_4.csv')

tobias_koop.head()
```

Note that statsmodels works seemlessly with Pandas dataframes (tobias_koop).

```{code-cell} ipython3
formula = "ln_wage ~ educ + pexp + pexp2 + broken_home"
results = smf.ols(formula,tobias_koop).fit()
print(results.summary())
```

Note that the results object contains the other types of information (including robust standard errors).  

```{code-cell} ipython3
results.HC0_se
```

If you want to specify which variance covariance matrix to use, use the `cov_type` argument.

```{code-cell} ipython3
results = smf.ols(formula,tobias_koop).fit(cov_type='HC0')
print(results.summary())
```

## Statistical Distributions

In the course, we will be evaluating cdf and pdf's, building likelihood and log-likelihood functions, taking random samples and perhaps other statistical operations.  We will primarily use the `scipy` family of distributions.  At times, I may "slip-up" and use some `numpy` functions, but in my opinion the `scipy` functions are more powerful by saving time coding.  Here are some examples:

### Evaluating pdf's and cdf's:


```{code-cell} ipython3
mean = 0
y = 1.5
std = 1
print("PDF value") 
print(norm(mean,std).pdf(y))
print("CDF value")
print(norm(mean,std).cdf(y))
```

Now plot the cumulative and probability density functions (CDF and PDF, respectively):

```{code-cell} ipython3
# use these for some plots
y_ = np.arange(-3.5,3.5,.05)

plot_pdf = [norm(mean, std).pdf(i) for i in y_]
plot_cdf = [norm(mean, std).cdf(i) for i in y_]

plt.figure(figsize=(8,8))
plt.plot(y_, plot_pdf, label="Standard Normal PDF")
plt.plot(y_, plot_cdf, label="Standard Normal CDF")
plt.legend(loc='upper left')
plt.ylim(0,1)
sbn.despine(offset=3.)
plt.show()
```

### Random Variates

Let's draw some random numbers from the normal distribution.  This time, instead of using mean 0 and standard deviation of 1, let's assume the distribution is centered at 10 with a standard deviation of 2.  Draw 100 random variates:


```{code-cell} ipython3
mean = 10
std = 2
N = 100
y = norm(mean,std).rvs(N)
print(y)
```

To convince ourselves, we have decent random numbers, let's view a histogram:

```{code-cell} ipython3
print(np.mean(y))
print(np.std(y))
plt.hist(y,bins=10)
plt.show()
```
 
As far as visually convincing us, maybe the jury is still out based on the histogram.  But with $N=100$, it is a small sample.

### Likelihood functions

Suppose the values y are observed and we want to calculate the likelihood of y given a mean for each of the 100 observations in y.  Even though we know the mean above should be 10, let's calculate the likelihood of y given a mean of 10.3 and a standard deviation of 2 assuming y is distributed normal.  Note these are simply pdf values for each datapoint:

$$
\mathcal{L}_i(y_i|\mu,\sigma) = \frac{1}{\sigma \sqrt{2\pi}}e^{\frac{(y_i - \mu)^2}{2\sigma^2}}
$$


```{code-cell} ipython3
norm(10.3,2).pdf(y)
```


We can also easily calculate the joint likelihood. This is

$$
\mathcal{L}(\mathbf{y}|\mu, \sigma) = \prod_{i \in N}\mathcal{L}_i(y_i|\mu,\sigma)
$$


```{code-cell} ipython3
norm(10.3,2).pdf(y).prod()
```

For computational advantages, we usually, work with log likelihoods, noting that

$$
ln(\mathcal{L}(\mathbf{y}|\mu, \sigma)) = \sum_{i \in N}ln(\mathcal{L}_i(y_i|\mu,\sigma))
$$

Implemented in code for our vector $\mathbf{y}$:


```{code-cell} ipython3
norm(10.3,2).logpdf(y).sum()
```

Let's visualize why the log-likelihood is better to work with than the likelihood:


```{code-cell} ipython3
mu_candidate = np.arange(7,13,.2) # plot likelihood and log-likelihood in this range
```


```{code-cell} ipython3
plot_pdf = np.array([norm(i,2.).pdf(y).prod() for i in mu_candidate])
plot_logpdf = np.array([norm(i,2.).logpdf(y).sum() for i in mu_candidate])

plt.figure(figsize=(10,5))
plt.subplot(121)
plt.title('Likelihood Function')
plt.plot(mu_candidate, plot_pdf, label="Likelihood Function")
plt.subplot(122)
plt.plot(mu_candidate, plot_logpdf, label="Log-Likelihood Function")
plt.title('Log-Likelihood Function')
sbn.despine()
plt.show()
```
    
## Maximum Likelihood Estimation

To find the maximum likelihood, solve the problem

   $$
   \max_{\mu, \sigma} \mathcal{L}(\mathbf{y}|\mu,\sigma) 
   $$

Oftentimes we express this using the log transformed version for the
reasons outlined above:
   
   $$
   \max_{\mu, \sigma} log(\mathcal{L}(\mathbf{y}|\mu,\sigma)) = \sum_{i=1}^N log(\mathcal{L}_i(y_i|\mu,\sigma))
   $$
This is termed the log-likelihood function.

To solve for estimates of $\mu$ and $\sigma$, either take derivatives and solve analytically, or use
numerical methods.
   
```{note}
The likelihood function tells us the **joint likelihood** of observing the sample $\mathbf{y}$ given estimates of $\mu$ and $\sigma$. 
```

This particular problem can (and probably should) be solved using analytical
methods, but since we will be using numerical methods primarily in
this course, we will demonstrate that in what follows.  

Looking at the numerical optimization routines available in `scipy`
you will note there are no entries for maximization, rather only
minimization (`scipy.optimize.minimize`).  We can use `minimize` for
finding the maximum of our log-likelihood function if we multiply it
by (-1).  So we define the negative log-likelihood as:

```{code-cell} ipython3
def neg_log_like(theta, y):
     # theta is our parameter vector
     #   parse into mu and sigma
     mu = theta[0]
     sigma = theta[-1]
     return -1.*norm(mu, sigma).logpdf(y).sum()
```

Then we can find the maximum likelihood estimate as:

```{code-cell} ipython3
startvals = np.array([10., 2.])
results = minimize(neg_log_like, startvals, args=(y))
results
```

Where the maximum likelihood estimates are 

```{code-cell} ipython3
print("MLE estimate for mu: ", results.x[0])
print("MLE estimate for sigma: ", results.x[1])
```

We can use the returned hessian inverse to compute standard errors.
We won't be doing any further maximum likelihood estimation in this
class.

```{note}
`scipy.minimize` can perform quite poorly with starting values that are far from the maximum likelihood values.  For these cases, use the `method='Nelder-Mead'` option.
```