---
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
---

# Introduction to Bayesian Econometrics

This course will cover the application of Bayesian statistical methods
for econometric inference.  Broadly speaking we will

1. Briefly discuss sampling methods for classical statistics
2. Introduce Bayes Rule and Provide an Application 
3. Examine the use of Monte Carlo Markov Chains
    * Link to Bayes Rule
    * Metropolis-Hastings and other Samplers
    * Chain "convergence" and diagnostics
4. Application: OLS, Time-Series Econometrics, Heirarchical Models

The **frequentist paradigm**, arguably the dominant statistical
paradigm in the social sciences (and what you all have studied),
relies on the following notions:

* $\beta$ is not random but neither is it known.  It is a fixed
  quantity.
* To uncover information about $\beta$, we observe part of some
  process (e.g. $\mathbf{y=x\beta+\epsilon}$).
* For statistical inference, we rely on **repeated trials** of
  $\mathbf{y}$ and $\mathbf{x}$, even if this repetition rarely (if
  ever) occurs in the social science context.
* $\mathbf{y}$ and $\mathbf{x}$ are considered random
* The model typically attempts to uncover information about
$\mathbf{\beta}$ by examining the likelihood function 
  
  $$
  prob(\mathbf{y}|\mathbf{b},\mathbf{x})
  $$ 
  
  where $\mathbf{b}$ are our estimates of $\beta$

The **bayesian paradigm** tackles the issue of estimating $\beta$ by

* Treating $\beta$ as random and unknown
* For a process ($\mathbf{y=x\beta+\epsilon}$) uncover information about $\beta$
* Treating $\mathbf{y}$ and $\mathbf{x}$ as fixed and non-random (at
  least once they are recorded in your dataset)
* Uncovers information about $\mathbf{\beta}$ by examining the
posterior likelihood 
  
  $$
  prob(\mathbf{b}|\mathbf{y},\mathbf{x})
  $$
  
  where $\mathbf{b}$ are our estimates of $\beta$

In quite a lot of instances, these two approaches give you the same
estimate for $\beta$.  Until recently, Bayesian Statistical modeling
wasn't used because calculating the posterior likelihood was
computationally challenging, but recent advances in the theory and
construction of Monte Carlo Markov Chains and computational ability
has really opened the door for Bayesian analysis for problems that
might not be estimated using the frequentist paradigm (ie. Maximim
Likelihood).  There is an **ongoing holy war** in the two statistical
camps, during the semester I will attempt to highlite the pros and
cons of each paradigm without taking a position on which one is
better.  My philosophy is that if it gets the job done, use it while
being aware of limitations and advantages.

## Repeated trials in frequentist statistics
A good jumping off point for this course is to understand the use of sampling techniques in a classical statistical paradigm.  
- Underlying all statistical inference that you have learned in statistics and econometrics is the idea of **repeated trials**.  Bootstrapping highlites this really well.  
- We will begin with an exploration of bootstrapping and see the
  implementation steps.

```{code-cell} ipython3
# load python libraries for this ipython notebook:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sbn
import statsmodels.formula.api as smf
import warnings
warnings.filterwarnings('ignore')

plt.style.use('ggplot')

np.random.seed(12578)
```

### Load Tobias and Koop data for `time==4`

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

And a data description,

```{code-cell} ipython3
tobias_koop.describe()
```

Next, visualize the distribution of `ln_wage` data, our "dependent
variable" for this problem.

```{code-cell} ipython3
plt.figure(figsize=(8,6))
tobias_koop.ln_wage.plot(kind='hist',bins=30)
plt.xlabel('ln_wage')
sbn.despine(offset=3.)
plt.show()
```

And using `statsmodels`, we can run this OLS regression:

$$
ln\_wage_i = \beta_0 + \beta_1 educ_i + \beta_2 pexp_i + \beta_3 broken\_home_i + \epsilon_i
$$

```{code-cell} ipython3
formula = 'ln_wage ~ educ + pexp + broken_home'

mod = smf.ols(formula, data=tobias_koop)
res=mod.fit()
res.summary()
```

These standard errors aren't robust (although we could ask for
that). As we saw in Cross-Section, we can non-parametrically bootstrap
confidence intervals of our parameters. While bootstrapping is
implemented in many ols routines, we implement them manually to see
how it works:

Steps for bootstrapping:

1. Loop $R$ times ($R$ = # of boostrapped replicates)
2. For each $r\in R$, sample $N$ times with replacement from rows in
   $\mathbf{Y}$ and $\mathbf{X}$, noting that $N$ is equal to the
   number of observations in the original dataset $\mathbf{Y}$ and
   $\mathbf{X}$. Denote each of these samples as $\mathbf{X}_r$ and
   $\mathbf{Y}_r$.
3. Estimate an OLS model and recover estimate for replicate $r$:
   $\hat{\beta}_r =
   (\mathbf{X}_r'\mathbf{X}_r)^{-1}\mathbf{X}_r'\mathbf{Y}_r$
4. Store results
5. Calculate the standard deviation (our estimate for the standard
   errors) of each parameter estimate, or construct non-parametric
   confidence intervals (preferred)

Note, for our purposes we don't need to see the regression results.
We only need to store them and move on to the next replicate.  Let's
investigate further the `res` object we created above.

### How to Recover OLS Estimates from Statsmodels

```{code-cell} ipython3
# The parameter values can be accessed here
beta = res.params
print(beta)
print("\nExtract the education parameter")
print(beta['educ'])
```

### Sampling with replacement

Fundamental to the idea of bootstrapping is to sample with replacement
to *preserve the shape of the underlying distribution of our
parameters*.  We can do this quite easily using numpy tools.

To do this, suppose we want to sample with replacement from 5 rows of
data indexed by `[0,1,2,3,4]`.

So we can sample with replacement from row numbers (index) and then
use that to select the rows we need for that replicate.  For each
replicate, run the model or calculate the statistic we are interested
in, store the results, and then repeat for all subsequent replicates.

```{code-cell} ipython3
row_id = np.arange(5)
print(row_id)

for r in range(3):
    this_sample = np.random.choice(row_id,size=row_id.shape[0],replace=True)
    print("\nReplicate", r+1)
    print("Rows to sample for this replicate")
    print(this_sample)
```

## The Bootstrap for our OLS Model

```{note}
Note that the boostrap is a method leaning on the frequentist idea of *sampling from the data* for inference about unknown parameters.  In this case we will be using it to learn about possible heterogeneity in individual-specific error standard variances $\sigma^2_i$ that is systematically related to our independent variables $\mathbf{x}$.
```

The following code uses the ideas above to perform the following steps
for each replicate $r$ of $R$ total replicates:
1. Sample with replacement from *rows* of our dataset
2. Calculate the statistic of interest (our OLS parameter estimates)
3. Store results

```{code-cell} ipython3
# Number of replicates
R = 2000 

# store each r of the R replicates here in rows:
results_boot = np.zeros((R,res.params.shape[0]))

row_id = range(0,tobias_koop.shape[0])

for r in range(R):
    # this samples with replacement from rows in the Tobias and Koop dataset
    this_sample = np.random.choice(row_id, size=tobias_koop.shape[0], 
                                   replace=True) # gives sampled row numbers
    # Define data for this replicate: 
    tobias_koop_r = tobias_koop.iloc[this_sample]
    # Estimate model
    results_r = smf.ols(formula,data=tobias_koop_r).fit().params
    # Store in row r of results_boot:
    results_boot[r,:] = results_r 
```

Let's see what is being stored in `results_boot`:
```{code-cell} ipython3
# Convert results to pandas dataframe for easier analysis:
results_boot = pd.DataFrame(results_boot,columns=['b_Intercept', 'b_educ',
                                                  'b_pexp', 'b_broken_home'])
results_boot.head(10)
```

and looking at the summary statistics:
```{code-cell} ipython3
results_boot.describe(percentiles=[.005,.025,.05,.5,.95,.975,.995])
```

From the above results, we see that the 95% confidence interval for
education is in the range [0.0668,0.10517] (approximately, as these
numbers will change each time you run the code above).  It is tempting
to interpret the confidence intervals above as **There is a 95% chance
that $\beta$ is in the range [0.0668,0.10517]**.  However, recall that
in the frequentist paradigm $\beta$ is fixed and non-random.  So it is
either in the range [0.0668,0.10517] or it isn't.  What is random is
the range [0.0668,0.10517].  So a better way to verbalize the
confidence interval is that if you repeated your regression analysis
many times, 95% of your calculated confidence interval would contain
the true parameter $\beta$.

Let's plot our bootstrap replicates, by parameter.  Note that the
values of each of these replicates are completely independent (or
should be).  Replicate 3 doesn't in any way inform us about
replicate 4.  For each parameter, we show the estimated parameter
(solid black line) and two parameter estimates:

1. The standard 95% confidence interval (dashed red line):

   $$
   CI_i = \hat{\beta}_i \pm 1.96 * \hat{\sigma}_{\beta_i} 
   $$

2. The 95% confidence interval calculated based off the 2.5 and 97.5
   percentiles from our boostrap database (dashed blue line).

### Focusing on `educ`

```{code-cell} ipython3
plt.figure(figsize=(10,8))
plt.xlabel('b_educ')
lw = 2
plt.axvline(beta['educ'],color='k',linestyle='dashed',lw=lw,
            label='OLS Estimate')
# Method 1. 95% Parametric CI's
plt.axvline(beta['educ'] -1.96*np.std(results_boot.b_educ),color='b',
            linestyle='dashed',
            lw=lw,label='Par. Lower 95%')
plt.axvline(beta['educ'] +1.96*np.std(results_boot.b_educ),color='b',
            linestyle='dashed',
            lw=lw,label='Par. Lower 95%')
# Method 2. Non-Parametric 95% CI's
plt.axvline(np.percentile(results_boot.b_educ,2.5),color='r',
            linestyle='dashed',
            lw=lw,label='Non-Par. Lower 95%')
plt.axvline(np.percentile(results_boot.b_educ,97.5),color='r',
            linestyle='dashed',
            lw=lw,label='Non-Par. Upper Lower 95%')
# scootch the upper limit of the x-axis a bit to the right for 
# non-overlapping legend
plt.xlim([.05,.15])

results_boot.b_educ.plot(kind='hist',bins=20,alpha=.5)
sbn.despine(offset=3.)
plt.legend()
plt.show()
```

We can see similar information for all parameters:

```{code-cell} ipython3
# set x axis values
replicate= np.arange(R)

# plot point estimate and confidence intervals:
plt.figure(figsize=(14, 20), dpi=200)
lw = 1

plt.subplot(4,2,1)
plt.ylabel("Intercept")
plt.plot(replicate, results_boot.b_Intercept, label="Intercept", lw=lw)
plt.subplot(4,2,2)
plt.hist(results_boot.b_Intercept,lw=lw, label="b_intercept", orientation='horizontal')

plt.subplot(4,2,3)
plt.ylabel("b_educ")
plt.plot(replicate, results_boot.b_educ, label="b_educ", lw=lw)
plt.subplot(4,2,4)
plt.hist(results_boot.b_educ,lw=lw, label="b_educ", orientation='horizontal')

plt.subplot(4,2,5)
plt.ylabel("b_pexp")
plt.plot(replicate, results_boot.b_pexp, label="b_pexp", lw=lw)
plt.subplot(4,2,6)
plt.hist(results_boot.b_pexp,lw=lw, label="b_pexp", orientation='horizontal')

plt.subplot(4,2,7)
plt.ylabel("b_broken_home")
plt.plot(replicate, results_boot.b_broken_home, label="b_broken_home", lw=lw)
plt.subplot(4,2,8)
plt.hist(results_boot.b_broken_home,lw=lw, label="b_broken_home", orientation='horizontal')
plt.show()
```

### Visualizing Covariances $\beta$:

```{code-cell} ipython3
plt.figure(figsize=(10,8));
sbn.jointplot(results_boot['b_educ'], results_boot['b_Intercept'], kind='hex');
plt.show()
```

or for all parameters:
```{code-cell} ipython3
g = sbn.PairGrid(results_boot)
g.map_diag(sbn.kdeplot)
g.map_offdiag(sbn.kdeplot, cmap="Blues_d")
plt.show()
```

```{note}
An important point about bootstrapping is that we are sampling with
replacement **from data** to uncover better information about our
error structure (something we need to estimate for the purposes of
inference).
```