Wooldridgeの練習問題をPythonで解く

# はじめに
- 最近、計量経済学の勉強を始めたので、教科書の練習問題をPythonで解きます。
- https://github.com/lystahi/wooldridge_python_notebook
- まだ途中なので、できていない部分や間違っている部分もあると思います。

# 教科書について
- _Introductory Econometrics: A Modern Approach_, by Jeffrey M. Wooldridge.
- https://www.amazon.co.jp/dp/1337558869/
- Rによる解法は[こちら](https://justinmshea.github.io/wooldridge/articles/Introductory-Econometrics-Examples.html)にあります。

# 解いていく

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
%matplotlib inline
```


```python
import warnings
warnings.simplefilter('ignore')
```

## Introduction

This notebook contains examples from _Introductory Econometrics: A Modern Approach, 6e_ by Jeffrey M. Wooldridge. Each example illustrates how to load data, build econometric models, and compute estimates with **Python**.

Now, install and load the `wooldridge` package and lets get started!


```python
# !pip install wooldridge
```


```python
from wooldridge import *
```

## Chapter 2: The Simple Regression Model

### **`Example 2.10:` A Log Wage Equation**

Load the `wage1` data and check out the documentation.


```python
dataWoo('wage1', description=True)
```

    name of dataset: wage1
    no of variables: 24
    no of observations: 526
    
    +----------+---------------------------------+
    | variable | label                           |
    +----------+---------------------------------+
    | wage     | average hourly earnings         |
    | educ     | years of education              |
    | exper    | years potential experience      |
    | tenure   | years with current employer     |
    | nonwhite | =1 if nonwhite                  |
    | female   | =1 if female                    |
    | married  | =1 if married                   |
    | numdep   | number of dependents            |
    | smsa     | =1 if live in SMSA              |
    | northcen | =1 if live in north central U.S |
    | south    | =1 if live in southern region   |
    | west     | =1 if live in western region    |
    | construc | =1 if work in construc. indus.  |
    | ndurman  | =1 if in nondur. manuf. indus.  |
    | trcommpu | =1 if in trans, commun, pub ut  |
    | trade    | =1 if in wholesale or retail    |
    | services | =1 if in services indus.        |
    | profserv | =1 if in prof. serv. indus.     |
    | profocc  | =1 if in profess. occupation    |
    | clerocc  | =1 if in clerical occupation    |
    | servocc  | =1 if in service occupation     |
    | lwage    | log(wage)                       |
    | expersq  | exper^2                         |
    | tenursq  | tenure^2                        |
    +----------+---------------------------------+
    
    These are data from the 1976 Current Population Survey, collected by
    Henry Farber when he and I were colleagues at MIT in 1988.


The documentation indicates these are data from the 1976 Current Population Survey, collected by Henry Farber when he and Wooldridge were colleagues at MIT in 1988.

**$wage$:** average hourly earnings

**$educ$:** years of education

First, make a scatter-plot of the two variables and look for possible patterns in the relationship between them.


```python
wage1 = dataWoo('wage1')
```


```python
plt.scatter(wage1.educ, wage1.wage)
plt.title("Wages vs. Education, 1976")
plt.xlabel("years of education")
plt.ylabel("Hourly wages")
plt.show()
```


![output_9_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/fca97dbc-38c7-9f4f-d7da-8ea25372bf37.png)


It appears that _**on average**_, more years of education, leads to higher wages.

The example in the text investigates what the _**percentage**_ change between wages and education might be. So, we must use the $log($`wage`$)$.

Build a linear model to estimate the relationship between the _log of wage_ (`lwage`) and _education_ (`educ`).

$$\widehat{log(wage)} = \beta_0 + \beta_1educ$$


```python
y = wage1.lwage
x = wage1.educ

model = sm.OLS(x, y)
results = model.fit()
```

Print the `summary` of the results.


```python
print(results.summary())
```

                                     OLS Regression Results                                
    =======================================================================================
    Dep. Variable:                   educ   R-squared (uncentered):                   0.916
    Model:                            OLS   Adj. R-squared (uncentered):              0.916
    Method:                 Least Squares   F-statistic:                              5717.
    Date:                Sat, 09 May 2020   Prob (F-statistic):                   2.18e-284
    Time:                        15:36:54   Log-Likelihood:                         -1438.9
    No. Observations:                 526   AIC:                                      2880.
    Df Residuals:                     525   BIC:                                      2884.
    Df Model:                           1                                                  
    Covariance Type:            nonrobust                                                  
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    lwage          7.2081      0.095     75.608      0.000       7.021       7.395
    ==============================================================================
    Omnibus:                       12.119   Durbin-Watson:                   1.699
    Prob(Omnibus):                  0.002   Jarque-Bera (JB):               13.454
    Skew:                          -0.295   Prob(JB):                      0.00120
    Kurtosis:                       3.515   Cond. No.                         1.00
    ==============================================================================
    
    Warnings:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.


Plot the $log($`wage`$)$ vs `educ`, adding a line representing the least squares
fit.


```python
x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values


def abline(slope, intercept):
    """Plot a line from slope and intercept"""
    axes = plt.gca()
    x_vals = np.array(axes.get_xlim())
    y_vals = intercept + slope * x_vals
    plt.plot(x_vals, y_vals, '--')


plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()
```


![output_15_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/588c23cd-97bb-a2ee-9eb0-9c3900722012.png)


## Chapter 3: Multiple Regression Analysis: Estimation

### **`Example 3.2:` Hourly Wage Equation**

Check the documentation for variable information


```python
dataWoo('wage1', description=True)
```

    name of dataset: wage1
    no of variables: 24
    no of observations: 526
    
    +----------+---------------------------------+
    | variable | label                           |
    +----------+---------------------------------+
    | wage     | average hourly earnings         |
    | educ     | years of education              |
    | exper    | years potential experience      |
    | tenure   | years with current employer     |
    | nonwhite | =1 if nonwhite                  |
    | female   | =1 if female                    |
    | married  | =1 if married                   |
    | numdep   | number of dependents            |
    | smsa     | =1 if live in SMSA              |
    | northcen | =1 if live in north central U.S |
    | south    | =1 if live in southern region   |
    | west     | =1 if live in western region    |
    | construc | =1 if work in construc. indus.  |
    | ndurman  | =1 if in nondur. manuf. indus.  |
    | trcommpu | =1 if in trans, commun, pub ut  |
    | trade    | =1 if in wholesale or retail    |
    | services | =1 if in services indus.        |
    | profserv | =1 if in prof. serv. indus.     |
    | profocc  | =1 if in profess. occupation    |
    | clerocc  | =1 if in clerical occupation    |
    | servocc  | =1 if in service occupation     |
    | lwage    | log(wage)                       |
    | expersq  | exper^2                         |
    | tenursq  | tenure^2                        |
    +----------+---------------------------------+
    
    These are data from the 1976 Current Population Survey, collected by
    Henry Farber when he and I were colleagues at MIT in 1988.


**$lwage$:** log of the average hourly earnings

**$educ$:** years of education

**$exper$:** years of potential experience

**$tenutre$:** years with current employer

Plot the variables against `lwage` and compare their distributions
and slope ($\beta$) of the simple regression lines.



```python
y = wage1.lwage
x = wage1.educ

model = sm.OLS(x, y)
results = model.fit()

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values

plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()
```


![output_19_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/7e47adef-491a-1bbf-1d56-b804803c58fb.png)




```python
y = wage1.lwage
x = wage1.exper

model = sm.OLS(x, y)
results = model.fit()

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values

plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()
```


![output_20_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/0245bc1f-3ba1-5468-b486-b91c978f0d9b.png)



```python
y = wage1.lwage
x = wage1.tenure

model = sm.OLS(x, y)
results = model.fit()

x_with_const = sm.add_constant(x)

model = sm.OLS(x_with_const, y)
results = model.fit()

slope, intercept = results.params[1].values/100, results.params[0].values

plt.scatter(x, y)
abline(slope, intercept)
plt.title("A Log Wage Equation")
plt.xlabel("years of education")
plt.ylabel("log of average hourly wages")
plt.show()
```


![output_21_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/84eb29d0-9df3-0baa-5b89-825b11562a85.png)



Estimate the model regressing _educ_, _exper_, and _tenure_ against _log(wage)_.

$$\widehat{log(wage)} = \beta_0 + \beta_1educ + \beta_3exper + \beta_4tenure$$


```python
results = smf.ols("lwage ~ educ + exper + tenure", data=wage1).fit()
results.params
```




    Intercept    0.284360
    educ         0.092029
    exper        0.004121
    tenure       0.022067
    dtype: float64




```python
coef_df = pd.DataFrame(results.params)
coef_df.drop(index=['Intercept'], inplace=True)
coef_df.rename(columns={0: "coef"}, inplace=True)
```


```python
coef_df
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>coef</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>educ</td>
      <td>0.092029</td>
    </tr>
    <tr>
      <td>exper</td>
      <td>0.004121</td>
    </tr>
    <tr>
      <td>tenure</td>
      <td>0.022067</td>
    </tr>
  </tbody>
</table>
</div>




```python
coef_df.sort_values(by='coef', ascending=False).plot.bar()
```




    <matplotlib.axes._subplots.AxesSubplot at 0x1c1b4c9390>




![output_26_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/4cfd0551-3bd3-4141-5c0c-7cc76841bee7.png)



## Chapter 4: Multiple Regression Analysis: Inference

### **`Example 4.1` Hourly Wage Equation**

Using the same model estimated in **`example: 3.2`**, examine and compare the standard errors associated with each coefficient. Like the textbook, these are contained in parenthesis next to each associated coefficient.


```python
results = smf.ols("lwage ~ educ + exper + tenure", data=wage1).fit()
results.summary()
```


<div>

<table class="simpletable1">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>          <td>lwage</td>      <th>  R-squared:         </th> <td>   0.316</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.312</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   80.39</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>9.13e-43</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:55</td>     <th>  Log-Likelihood:    </th> <td> -313.55</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   526</td>      <th>  AIC:               </th> <td>   635.1</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   522</td>      <th>  BIC:               </th> <td>   652.2</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable2">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    0.2844</td> <td>    0.104</td> <td>    2.729</td> <td> 0.007</td> <td>    0.080</td> <td>    0.489</td>
</tr>
<tr>
  <th>educ</th>      <td>    0.0920</td> <td>    0.007</td> <td>   12.555</td> <td> 0.000</td> <td>    0.078</td> <td>    0.106</td>
</tr>
<tr>
  <th>exper</th>     <td>    0.0041</td> <td>    0.002</td> <td>    2.391</td> <td> 0.017</td> <td>    0.001</td> <td>    0.008</td>
</tr>
<tr>
  <th>tenure</th>    <td>    0.0221</td> <td>    0.003</td> <td>    7.133</td> <td> 0.000</td> <td>    0.016</td> <td>    0.028</td>
</tr>
</table>
<table class="simpletable3">
<tr>
  <th>Omnibus:</th>       <td>11.534</td> <th>  Durbin-Watson:     </th> <td>   1.769</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.003</td> <th>  Jarque-Bera (JB):  </th> <td>  20.941</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 0.021</td> <th>  Prob(JB):          </th> <td>2.84e-05</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.977</td> <th>  Cond. No.          </th> <td>    135.</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>


```python
results.params[2]/results.bse[2]
```




    2.3914371073142298



### **`Example 4.7` Effect of Job Training on Firm Scrap Rates**

Load the `jtrain` data set. 


```python
dataWoo('jtrain', description=True)
```

    name of dataset: jtrain
    no of variables: 30
    no of observations: 471
    
    +----------+---------------------------------+
    | variable | label                           |
    +----------+---------------------------------+
    | year     | 1987, 1988, or 1989             |
    | fcode    | firm code number                |
    | employ   | # employees at plant            |
    | sales    | annual sales, $                 |
    | avgsal   | average employee salary         |
    | scrap    | scrap rate (per 100 items)      |
    | rework   | rework rate (per 100 items)     |
    | tothrs   | total hours training            |
    | union    | =1 if unionized                 |
    | grant    | = 1 if received grant           |
    | d89      | = 1 if year = 1989              |
    | d88      | = 1 if year = 1988              |
    | totrain  | total employees trained         |
    | hrsemp   | tothrs/totrain                  |
    | lscrap   | log(scrap)                      |
    | lemploy  | log(employ)                     |
    | lsales   | log(sales)                      |
    | lrework  | log(rework)                     |
    | lhrsemp  | log(1 + hrsemp)                 |
    | lscrap_1 | lagged lscrap; missing 1987     |
    | grant_1  | lagged grant; assumed 0 in 1987 |
    | clscrap  | lscrap - lscrap_1; year > 1987  |
    | cgrant   | grant - grant_1                 |
    | clemploy | lemploy - lemploy[_n-1]         |
    | clsales  | lavgsal - lavgsal[_n-1]         |
    | lavgsal  | log(avgsal)                     |
    | clavgsal | lavgsal - lavgsal[_n-1]         |
    | cgrant_1 | cgrant[_n-1]                    |
    | chrsemp  | hrsemp - hrsemp[_n-1]           |
    | clhrsemp | lhrsemp - lhrsemp[_n-1]         |
    +----------+---------------------------------+
    
    H. Holzer, R. Block, M. Cheatham, and J. Knott (1993), “Are Training
    Subsidies Effective? The Michigan Experience,” Industrial and Labor
    Relations Review 46, 625-636. The authors kindly provided the data.



```python
jtrain = dataWoo('jtrain')
jtrain_subset = jtrain[jtrain.year == 1987][['year', 'union', 'lscrap',
                                             'hrsemp', 'lsales', 'lemploy']]
jtrain_subset.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thread>
    <tr style="text-align: right;">
      <th></th>
      <th>year</th>
      <th>union</th>
      <th>lscrap</th>
      <th>hrsemp</th>
      <th>lsales</th>
      <th>lemploy</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>0</td>
      <td>1987</td>
      <td>0</td>
      <td>NaN</td>
      <td>12.0</td>
      <td>17.665659</td>
      <td>4.605170</td>
    </tr>
    <tr>
      <td>3</td>
      <td>1987</td>
      <td>0</td>
      <td>NaN</td>
      <td>12.0</td>
      <td>14.260197</td>
      <td>2.484907</td>
    </tr>
    <tr>
      <td>6</td>
      <td>1987</td>
      <td>0</td>
      <td>NaN</td>
      <td>37.5</td>
      <td>13.527828</td>
      <td>2.995732</td>
    </tr>
    <tr>
      <td>9</td>
      <td>1987</td>
      <td>0</td>
      <td>NaN</td>
      <td>0.0</td>
      <td>16.982714</td>
      <td>5.298317</td>
    </tr>
    <tr>
      <td>12</td>
      <td>1987</td>
      <td>0</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>15.607270</td>
      <td>NaN</td>
    </tr>
  </tbody>
</table>
</div>




```python
jtrain_subset.isnull().sum()
```




    year         0
    union        0
    lscrap     103
    hrsemp      28
    lsales      38
    lemploy     13
    dtype: int64




```python
jtrain_clean = jtrain_subset.dropna()
```

Now create the linear model regressing `hrsemp`(total hours training/total employees trained), `lsales`(log of annual sales), and `lemploy`(the log of the number of the employees), against `lscrap`(the log of the scrape rate).

$$lscrap = \alpha + \beta_1 hrsemp + \beta_2 lsales + \beta_3 lemploy$$


```python
results = smf.ols("lscrap ~ hrsemp + lsales + lemploy",
                  data=jtrain_clean).fit()
results.summary()
```


<div>

<table class="simpletable4">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lscrap</td>      <th>  R-squared:         </th> <td>   0.310</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.257</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   5.838</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th>  <td>0.00215</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:55</td>     <th>  Log-Likelihood:    </th> <td> -70.198</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>    43</td>      <th>  AIC:               </th> <td>   148.4</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>    39</td>      <th>  BIC:               </th> <td>   155.4</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable5">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>   11.7443</td> <td>    4.575</td> <td>    2.567</td> <td> 0.014</td> <td>    2.491</td> <td>   20.997</td>
</tr>
<tr>
  <th>hrsemp</th>    <td>   -0.0422</td> <td>    0.019</td> <td>   -2.259</td> <td> 0.030</td> <td>   -0.080</td> <td>   -0.004</td>
</tr>
<tr>
  <th>lsales</th>    <td>   -0.9506</td> <td>    0.370</td> <td>   -2.570</td> <td> 0.014</td> <td>   -1.699</td> <td>   -0.203</td>
</tr>
<tr>
  <th>lemploy</th>   <td>    0.9921</td> <td>    0.357</td> <td>    2.780</td> <td> 0.008</td> <td>    0.270</td> <td>    1.714</td>
</tr>
</table>
<table class="simpletable6">
<tr>
  <th>Omnibus:</th>       <td> 1.109</td> <th>  Durbin-Watson:     </th> <td>   2.163</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.574</td> <th>  Jarque-Bera (JB):  </th> <td>   0.934</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 0.093</td> <th>  Prob(JB):          </th> <td>   0.627</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 2.302</td> <th>  Cond. No.          </th> <td>    417.</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>

## Chapter 5: Multiple Regression Analysis: OLS Asymptotics

### **`Example 5.1:` Housing Prices and Distance From an Incinerator**

Load the `hprice3` data set.


```python
dataWoo('hprice3', description=True)
```

    name of dataset: hprice3
    no of variables: 19
    no of observations: 321
    
    +----------+---------------------------------+
    | variable | label                           |
    +----------+---------------------------------+
    | year     | 1978, 1981                      |
    | age      | age of house                    |
    | agesq    | age^2                           |
    | nbh      | neighborhood, 1-6               |
    | cbd      | dist. to cent. bus. dstrct, ft. |
    | inst     | dist. to interstate, ft.        |
    | linst    | log(inst)                       |
    | price    | selling price                   |
    | rooms    | # rooms in house                |
    | area     | square footage of house         |
    | land     | square footage lot              |
    | baths    | # bathrooms                     |
    | dist     | dist. from house to incin., ft. |
    | ldist    | log(dist)                       |
    | lprice   | log(price)                      |
    | y81      | =1 if year = 1981               |
    | larea    | log(area)                       |
    | lland    | log(land)                       |
    | linstsq  | linst^2                         |
    +----------+---------------------------------+
    
    -



```python
hprice3 = dataWoo('hprice3')
```

Next, model the $log($`price`$)$ against the $log($`dist`$)$ to estimate the percentage relationship between the two.


```python
price_dist_model = smf.ols("lprice ~ ldist", data=hprice3).fit()
```

Create another model that controls for "quality" variables, such as square footage `area` per house.


```python
price_area_model = smf.ols("lprice ~ ldist + larea", data=hprice3).fit()
```


```python
price_dist_model.summary()
```



<div>
<table class="simpletable7">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lprice</td>      <th>  R-squared:         </th> <td>   0.120</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.117</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   43.48</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>1.78e-10</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:55</td>     <th>  Log-Likelihood:    </th> <td> -169.60</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   321</td>      <th>  AIC:               </th> <td>   343.2</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   319</td>      <th>  BIC:               </th> <td>   350.7</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable8">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    8.2575</td> <td>    0.474</td> <td>   17.427</td> <td> 0.000</td> <td>    7.325</td> <td>    9.190</td>
</tr>
<tr>
  <th>ldist</th>     <td>    0.3172</td> <td>    0.048</td> <td>    6.594</td> <td> 0.000</td> <td>    0.223</td> <td>    0.412</td>
</tr>
</table>
<table class="simpletable9">
<tr>
  <th>Omnibus:</th>       <td> 3.073</td> <th>  Durbin-Watson:     </th> <td>   0.944</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.215</td> <th>  Jarque-Bera (JB):  </th> <td>   2.868</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 0.228</td> <th>  Prob(JB):          </th> <td>   0.238</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.079</td> <th>  Cond. No.          </th> <td>    205.</td>
</tr>
</table><br/><br/>Warnings:<br/><p>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.</p>

</div>




```python
price_area_model.summary()
```



<div>
<table class="simpletable10">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lprice</td>      <th>  R-squared:         </th> <td>   0.474</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.471</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   143.2</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>4.53e-45</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:55</td>     <th>  Log-Likelihood:    </th> <td> -87.041</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   321</td>      <th>  AIC:               </th> <td>   180.1</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   318</td>      <th>  BIC:               </th> <td>   191.4</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable11">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    3.4939</td> <td>    0.491</td> <td>    7.121</td> <td> 0.000</td> <td>    2.529</td> <td>    4.459</td>
</tr>
<tr>
  <th>ldist</th>     <td>    0.1962</td> <td>    0.038</td> <td>    5.142</td> <td> 0.000</td> <td>    0.121</td> <td>    0.271</td>
</tr>
<tr>
  <th>larea</th>     <td>    0.7837</td> <td>    0.054</td> <td>   14.625</td> <td> 0.000</td> <td>    0.678</td> <td>    0.889</td>
</tr>
</table>
<table class="simpletable12">
<tr>
  <th>Omnibus:</th>       <td>13.639</td> <th>  Durbin-Watson:     </th> <td>   1.118</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.001</td> <th>  Jarque-Bera (JB):  </th> <td>  18.845</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.334</td> <th>  Prob(JB):          </th> <td>8.09e-05</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.982</td> <th>  Cond. No.          </th> <td>    345.</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
</div>


## Chapter 6: Multiple Regression: Further Issues

### **`Example 6.1:` Effects of Pollution on Housing Prices, standardized.**

Load the `hprice2` data and view the documentation.


```python
dataWoo('hprice2', description=True)
```

    name of dataset: hprice2
    no of variables: 12
    no of observations: 506
    
    +----------+-------------------------------+
    | variable | label                         |
    +----------+-------------------------------+
    | price    | median housing price, $       |
    | crime    | crimes committed per capita   |
    | nox      | nit ox concen; parts per 100m |
    | rooms    | avg number of rooms           |
    | dist     | wght dist to 5 employ centers |
    | radial   | access. index to rad. hghwys  |
    | proptax  | property tax per $1000        |
    | stratio  | average student-teacher ratio |
    | lowstat  | perc of people 'lower status' |
    | lprice   | log(price)                    |
    | lnox     | log(nox)                      |
    | lproptax | log(proptax)                  |
    +----------+-------------------------------+
    
    D. Harrison and D.L. Rubinfeld (1978), “Hedonic Housing Prices and the
    Demand for Clean Air,” by Harrison, D. and D.L.Rubinfeld, Journal of
    Environmental Economics and Management 5, 81-102. Diego Garcia, a
    former Ph.D. student in economics at MIT, kindly provided these data,
    which he obtained from the book Regression Diagnostics: Identifying
    Influential Data and Sources of Collinearity, by D.A. Belsey, E. Kuh,
    and R. Welsch, 1990. New York: Wiley.


$$price = \beta_0 + \beta_1nox + \beta_2crime + \beta_3rooms + \beta_4dist + \beta_5stratio + \mu$$

Estimate the usual `lm` model.


```python
hprice2 = dataWoo('hprice2')
housing_level = smf.ols("price ~ nox + crime + rooms + dist + stratio",
                        data=hprice2).fit()
```

Estimate the same model, but standardized coefficients by wrapping each variable with R’s scale function:

$$\widehat{zprice} = \beta_1znox + \beta_2zcrime + \beta_3zrooms + \beta_4zdist + \beta_5zstratio$$


```python
from sklearn.preprocessing import scale

zprice = pd.DataFrame({"zprice":scale(hprice2["price"])})
znox = pd.DataFrame({"znox":scale(hprice2["nox"])})
zcrime = pd.DataFrame({"zcrime":scale(hprice2["crime"])})
zrooms = pd.DataFrame({"zrooms":scale(hprice2["rooms"])})
zdist = pd.DataFrame({"zdist":scale(hprice2["dist"])})
zstratio = pd.DataFrame({"zstratio":scale(hprice2["stratio"])})
hprice2 = pd.concat([hprice2, zprice, znox, zcrime, zrooms, zdist, zstratio],axis=1)

housing_standardized = smf.ols("zprice ~ znox + zcrime + zrooms + zdist + zstratio"
                               , data = hprice2).fit()
```


```python
housing_level.summary()
```



<div>
<table class="simpletable13">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>          <td>price</td>      <th>  R-squared:         </th> <td>   0.636</td> 
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.632</td> 
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   174.5</td> 
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>3.61e-107</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -5080.8</td> 
</tr>
<tr>
  <th>No. Observations:</th>      <td>   506</td>      <th>  AIC:               </th> <td>1.017e+04</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   500</td>      <th>  BIC:               </th> <td>1.020e+04</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     5</td>      <th>                     </th>     <td> </td>    
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    
</tr>
</table>
<table class="simpletable14">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td> 2.087e+04</td> <td> 5054.599</td> <td>    4.129</td> <td> 0.000</td> <td> 1.09e+04</td> <td> 3.08e+04</td>
</tr>
<tr>
  <th>nox</th>       <td>-2706.4326</td> <td>  354.087</td> <td>   -7.643</td> <td> 0.000</td> <td>-3402.114</td> <td>-2010.751</td>
</tr>
<tr>
  <th>crime</th>     <td> -153.6010</td> <td>   32.929</td> <td>   -4.665</td> <td> 0.000</td> <td> -218.297</td> <td>  -88.905</td>
</tr>
<tr>
  <th>rooms</th>     <td> 6735.4983</td> <td>  393.604</td> <td>   17.112</td> <td> 0.000</td> <td> 5962.177</td> <td> 7508.819</td>
</tr>
<tr>
  <th>dist</th>      <td>-1026.8063</td> <td>  188.108</td> <td>   -5.459</td> <td> 0.000</td> <td>-1396.386</td> <td> -657.227</td>
</tr>
<tr>
  <th>stratio</th>   <td>-1149.2038</td> <td>  127.429</td> <td>   -9.018</td> <td> 0.000</td> <td>-1399.566</td> <td> -898.842</td>
</tr>
</table>
<table class="simpletable15">
<tr>
  <th>Omnibus:</th>       <td>272.145</td> <th>  Durbin-Watson:     </th> <td>   0.865</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>2647.578</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 2.150</td>  <th>  Prob(JB):          </th> <td>    0.00</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td>13.348</td>  <th>  Cond. No.          </th> <td>    432.</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
</div>



```python
housing_standardized.summary()
```



<div>
<table class="simpletable16">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>zprice</td>      <th>  R-squared:         </th> <td>   0.636</td> 
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.632</td> 
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   174.5</td> 
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>3.61e-107</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -462.53</td> 
</tr>
<tr>
  <th>No. Observations:</th>      <td>   506</td>      <th>  AIC:               </th> <td>   937.1</td> 
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   500</td>      <th>  BIC:               </th> <td>   962.4</td> 
</tr>
<tr>
  <th>Df Model:</th>              <td>     5</td>      <th>                     </th>     <td> </td>    
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    
</tr>
</table>
<table class="simpletable17">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td> 9.368e-17</td> <td>    0.027</td> <td> 3.47e-15</td> <td> 1.000</td> <td>   -0.053</td> <td>    0.053</td>
</tr>
<tr>
  <th>znox</th>      <td>   -0.3404</td> <td>    0.045</td> <td>   -7.643</td> <td> 0.000</td> <td>   -0.428</td> <td>   -0.253</td>
</tr>
<tr>
  <th>zcrime</th>    <td>   -0.1433</td> <td>    0.031</td> <td>   -4.665</td> <td> 0.000</td> <td>   -0.204</td> <td>   -0.083</td>
</tr>
<tr>
  <th>zrooms</th>    <td>    0.5139</td> <td>    0.030</td> <td>   17.112</td> <td> 0.000</td> <td>    0.455</td> <td>    0.573</td>
</tr>
<tr>
  <th>zdist</th>     <td>   -0.2348</td> <td>    0.043</td> <td>   -5.459</td> <td> 0.000</td> <td>   -0.319</td> <td>   -0.150</td>
</tr>
<tr>
  <th>zstratio</th>  <td>   -0.2703</td> <td>    0.030</td> <td>   -9.018</td> <td> 0.000</td> <td>   -0.329</td> <td>   -0.211</td>
</tr>
</table>
<table class="simpletable18">
<tr>
  <th>Omnibus:</th>       <td>272.145</td> <th>  Durbin-Watson:     </th> <td>   0.865</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>2647.578</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 2.150</td>  <th>  Prob(JB):          </th> <td>    0.00</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td>13.348</td>  <th>  Cond. No.          </th> <td>    3.33</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>

### **`Example 6.2:` Effects of Pollution on Housing Prices, Quadratic Interactive Term**

Modify the housing model from **`example 4.5`**, adding a quadratic term in _rooms_: 

$$log(price) = \beta_0 + \beta_1log(nox) + \beta_2log(dist) + \beta_3rooms + \beta_4rooms^2 + \beta_5stratio + \mu$$


```python
ldist = pd.DataFrame({"ldist":np.log(hprice2["dist"])})
roomssq = pd.DataFrame({"roomssq":hprice2["rooms"]**2})
hprice2 = pd.concat([hprice2, ldist, roomssq],axis=1)

housing_model_4_5 = smf.ols("lprice ~ lnox + ldist + rooms + stratio"
                            , data=hprice2).fit()
housing_model_6_2 = smf.ols("lprice ~ lnox + ldist + rooms + roomssq + stratio"
                            , data=hprice2).fit()
```


```python
housing_model_4_5.summary()
```


<div>

<table class="simpletable19">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lprice</td>      <th>  R-squared:         </th> <td>   0.584</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.581</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   175.9</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>5.53e-94</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -43.495</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   506</td>      <th>  AIC:               </th> <td>   96.99</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   501</td>      <th>  BIC:               </th> <td>   118.1</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     4</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable20">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>   11.0839</td> <td>    0.318</td> <td>   34.843</td> <td> 0.000</td> <td>   10.459</td> <td>   11.709</td>
</tr>
<tr>
  <th>lnox</th>      <td>   -0.9535</td> <td>    0.117</td> <td>   -8.168</td> <td> 0.000</td> <td>   -1.183</td> <td>   -0.724</td>
</tr>
<tr>
  <th>ldist</th>     <td>   -0.1343</td> <td>    0.043</td> <td>   -3.117</td> <td> 0.002</td> <td>   -0.219</td> <td>   -0.050</td>
</tr>
<tr>
  <th>rooms</th>     <td>    0.2545</td> <td>    0.019</td> <td>   13.736</td> <td> 0.000</td> <td>    0.218</td> <td>    0.291</td>
</tr>
<tr>
  <th>stratio</th>   <td>   -0.0525</td> <td>    0.006</td> <td>   -8.894</td> <td> 0.000</td> <td>   -0.064</td> <td>   -0.041</td>
</tr>
</table>
<table class="simpletable21">
<tr>
  <th>Omnibus:</th>       <td>61.317</td> <th>  Durbin-Watson:     </th> <td>   0.682</td> 
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td> 480.143</td> 
</tr>
<tr>
  <th>Skew:</th>          <td> 0.051</td> <th>  Prob(JB):          </th> <td>5.47e-105</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 7.771</td> <th>  Cond. No.          </th> <td>    560.</td> 
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>


```python
housing_model_6_2.summary()
```


<div>

<table class="simpletable22">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lprice</td>      <th>  R-squared:         </th> <td>   0.603</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.599</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   151.8</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>7.89e-98</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -31.806</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   506</td>      <th>  AIC:               </th> <td>   75.61</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   500</td>      <th>  BIC:               </th> <td>   101.0</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     5</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>   13.3855</td> <td>    0.566</td> <td>   23.630</td> <td> 0.000</td> <td>   12.273</td> <td>   14.498</td>
</tr>
<tr>
  <th>lnox</th>      <td>   -0.9017</td> <td>    0.115</td> <td>   -7.862</td> <td> 0.000</td> <td>   -1.127</td> <td>   -0.676</td>
</tr>
<tr>
  <th>ldist</th>     <td>   -0.0868</td> <td>    0.043</td> <td>   -2.005</td> <td> 0.045</td> <td>   -0.172</td> <td>   -0.002</td>
</tr>
<tr>
  <th>rooms</th>     <td>   -0.5451</td> <td>    0.165</td> <td>   -3.295</td> <td> 0.001</td> <td>   -0.870</td> <td>   -0.220</td>
</tr>
<tr>
  <th>roomssq</th>   <td>    0.0623</td> <td>    0.013</td> <td>    4.862</td> <td> 0.000</td> <td>    0.037</td> <td>    0.087</td>
</tr>
<tr>
  <th>stratio</th>   <td>   -0.0476</td> <td>    0.006</td> <td>   -8.129</td> <td> 0.000</td> <td>   -0.059</td> <td>   -0.036</td>
</tr>
</table>
<table class="simpletable23">
<tr>
  <th>Omnibus:</th>       <td>56.649</td> <th>  Durbin-Watson:     </th> <td>   0.691</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td> 384.168</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.100</td> <th>  Prob(JB):          </th> <td>3.79e-84</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 7.264</td> <th>  Cond. No.          </th> <td>2.30e+03</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 2.3e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems.

</div>

## Chapter 7: Multiple Regression Analysis with Qualitative Information 

### **`Example 7.4:` Housing Price Regression, Qualitative Binary variable**

This time, use the `hrprice1` data. 


```python
dataWoo('hprice1', description=True)
```

    name of dataset: hprice1
    no of variables: 10
    no of observations: 88
    
    +----------+------------------------------+
    | variable | label                        |
    +----------+------------------------------+
    | price    | house price, $1000s          |
    | assess   | assessed value, $1000s       |
    | bdrms    | number of bdrms              |
    | lotsize  | size of lot in square feet   |
    | sqrft    | size of house in square feet |
    | colonial | =1 if home is colonial style |
    | lprice   | log(price)                   |
    | lassess  | log(assess                   |
    | llotsize | log(lotsize)                 |
    | lsqrft   | log(sqrft)                   |
    +----------+------------------------------+
    
    Collected from the real estate pages of the Boston Globe during 1990.
    These are homes that sold in the Boston, MA area.


$$\widehat{log(price)} = \beta_0 + \beta_1log(lotsize) + \beta_2log(sqrft) + \beta_3bdrms + \beta_4colonial $$

Estimate the coefficients of the above linear model on the `hprice` data set.


```python
hprice1 = dataWoo('hprice1')
housing_qualitative = smf.ols("lprice ~ llotsize + lsqrft + bdrms + colonial",
                              data=hprice1).fit()
housing_qualitative.summary()
```


<div>

<table class="simpletable24">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lprice</td>      <th>  R-squared:         </th> <td>   0.649</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.632</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   38.38</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>3.74e-18</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td>  26.619</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>    88</td>      <th>  AIC:               </th> <td>  -43.24</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>    83</td>      <th>  BIC:               </th> <td>  -30.85</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     4</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable25">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>   -1.3496</td> <td>    0.651</td> <td>   -2.073</td> <td> 0.041</td> <td>   -2.644</td> <td>   -0.055</td>
</tr>
<tr>
  <th>llotsize</th>  <td>    0.1678</td> <td>    0.038</td> <td>    4.395</td> <td> 0.000</td> <td>    0.092</td> <td>    0.244</td>
</tr>
<tr>
  <th>lsqrft</th>    <td>    0.7072</td> <td>    0.093</td> <td>    7.620</td> <td> 0.000</td> <td>    0.523</td> <td>    0.892</td>
</tr>
<tr>
  <th>bdrms</th>     <td>    0.0268</td> <td>    0.029</td> <td>    0.934</td> <td> 0.353</td> <td>   -0.030</td> <td>    0.084</td>
</tr>
<tr>
  <th>colonial</th>  <td>    0.0538</td> <td>    0.045</td> <td>    1.202</td> <td> 0.233</td> <td>   -0.035</td> <td>    0.143</td>
</tr>
</table>
<table class="simpletable26">
<tr>
  <th>Omnibus:</th>       <td>13.728</td> <th>  Durbin-Watson:     </th> <td>   2.077</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.001</td> <th>  Jarque-Bera (JB):  </th> <td>  50.828</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.053</td> <th>  Prob(JB):          </th> <td>9.18e-12</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 6.722</td> <th>  Cond. No.          </th> <td>    411.</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>

## Chapter 8: Heteroskedasticity

### **`Example 8.9:` Determinants of Personal Computer Ownership**

$$\widehat{PC} = \beta_0 + \beta_1hsGPA + \beta_2ACT + \beta_3parcoll + \beta_4colonial $$
Christopher Lemmon, a former MSU undergraduate, collected these data from a survey he took of MSU students in Fall 1994. Load `gpa1` and create a new variable combining the `fathcoll` and `mothcoll`, into `parcoll`. This new column indicates if either parent went to college.


```python
dataWoo('gpa1', description=True)
```

    name of dataset: gpa1
    no of variables: 29
    no of observations: 141
    
    +----------+--------------------------------+
    | variable | label                          |
    +----------+--------------------------------+
    | age      | in years                       |
    | soph     | =1 if sophomore                |
    | junior   | =1 if junior                   |
    | senior   | =1 if senior                   |
    | senior5  | =1 if fifth year senior        |
    | male     | =1 if male                     |
    | campus   | =1 if live on campus           |
    | business | =1 if business major           |
    | engineer | =1 if engineering major        |
    | colGPA   | MSU GPA                        |
    | hsGPA    | high school GPA                |
    | ACT      | 'achievement' score            |
    | job19    | =1 if job <= 19 hours          |
    | job20    | =1 if job >= 20 hours          |
    | drive    | =1 if drive to campus          |
    | bike     | =1 if bicycle to campus        |
    | walk     | =1 if walk to campus           |
    | voluntr  | =1 if do volunteer work        |
    | PC       | =1 of pers computer at sch     |
    | greek    | =1 if fraternity or sorority   |
    | car      | =1 if own car                  |
    | siblings | =1 if have siblings            |
    | bgfriend | =1 if boy- or girlfriend       |
    | clubs    | =1 if belong to MSU club       |
    | skipped  | avg lectures missed per week   |
    | alcohol  | avg # days per week drink alc. |
    | gradMI   | =1 if Michigan high school     |
    | fathcoll | =1 if father college grad      |
    | mothcoll | =1 if mother college grad      |
    +----------+--------------------------------+
    
    Christopher Lemmon, a former MSU undergraduate, collected these data
    from a survey he took of MSU students in Fall 1994.



```python
gpa1 = dataWoo('gpa1')

gpa1['parcoll'] = 0
gpa1.parcoll.loc[gpa1.fathcoll == 1 | gpa1.mothcoll] = 1
```


```python
GPA_OLS = smf.ols("PC ~ hsGPA + ACT + parcoll", data=gpa1).fit()
GPA_OLS.summary()
```


<div>

<table class="simpletable27">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>           <td>PC</td>        <th>  R-squared:         </th> <td>   0.023</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.002</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   1.087</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th>  <td> 0.357</td> 
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -97.630</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   141</td>      <th>  AIC:               </th> <td>   203.3</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   137</td>      <th>  BIC:               </th> <td>   215.1</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable28">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    0.1107</td> <td>    0.492</td> <td>    0.225</td> <td> 0.822</td> <td>   -0.862</td> <td>    1.083</td>
</tr>
<tr>
  <th>hsGPA</th>     <td>    0.0489</td> <td>    0.138</td> <td>    0.354</td> <td> 0.724</td> <td>   -0.224</td> <td>    0.322</td>
</tr>
<tr>
  <th>ACT</th>       <td>    0.0014</td> <td>    0.016</td> <td>    0.090</td> <td> 0.929</td> <td>   -0.030</td> <td>    0.032</td>
</tr>
<tr>
  <th>parcoll</th>   <td>    0.1463</td> <td>    0.085</td> <td>    1.728</td> <td> 0.086</td> <td>   -0.021</td> <td>    0.314</td>
</tr>
</table>
<table class="simpletable29">
<tr>
  <th>Omnibus:</th>       <td>1290.376</td> <th>  Durbin-Watson:     </th> <td>   1.780</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>  21.607</td>
</tr>
<tr>
  <th>Skew:</th>           <td> 0.405</td>  <th>  Prob(JB):          </th> <td>2.03e-05</td>
</tr>
<tr>
  <th>Kurtosis:</th>       <td> 1.262</td>  <th>  Cond. No.          </th> <td>    298.</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>

## Chapter 9: More on Specification and Data Issues

### **`Example 9.8:` R&D Intensity and Firm Size**


$$rdintens = \beta_0 + \beta_1sales + \beta_2profmarg + \mu$$

From _Businessweek R&D Scoreboard_, October 25, 1991. Load the data and estimate the model.


```python
rdchem = dataWoo("rdchem")

all_rdchem = smf.ols("rdintens ~ sales + profmarg", data=rdchem).fit()
```


```python
all_rdchem.summary()
```


<div>

<table class="simpletable30">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>        <td>rdintens</td>     <th>  R-squared:         </th> <td>   0.076</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.012</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   1.195</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th>  <td> 0.317</td> 
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -63.725</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>    32</td>      <th>  AIC:               </th> <td>   133.4</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>    29</td>      <th>  BIC:               </th> <td>   137.8</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable31">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    2.6253</td> <td>    0.586</td> <td>    4.484</td> <td> 0.000</td> <td>    1.428</td> <td>    3.823</td>
</tr>
<tr>
  <th>sales</th>     <td> 5.338e-05</td> <td> 4.41e-05</td> <td>    1.211</td> <td> 0.236</td> <td>-3.68e-05</td> <td>    0.000</td>
</tr>
<tr>
  <th>profmarg</th>  <td>    0.0446</td> <td>    0.046</td> <td>    0.966</td> <td> 0.342</td> <td>   -0.050</td> <td>    0.139</td>
</tr>
</table>
<table class="simpletable32">
<tr>
  <th>Omnibus:</th>       <td>20.499</td> <th>  Durbin-Watson:     </th> <td>   1.694</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  28.349</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 1.674</td> <th>  Prob(JB):          </th> <td>6.98e-07</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 6.170</td> <th>  Cond. No.          </th> <td>1.49e+04</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.49e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems.

</div>

So, we can estimate the model without that data point to gain a better understanding of how `sales` and `profmarg` describe `rdintens` for most firms. We can use the `subset` argument of the linear model function to indicate that we only want to estimate the model using data that is less than the highest sales.


```python
rdchem_sub = rdchem[rdchem.sales < max(rdchem.sales)]

smallest_rdchem = smf.ols("rdintens ~ sales + profmarg", data=rdchem_sub).fit()
```


```python
smallest_rdchem.summary()
```


<div>

<table class="simpletable33">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>        <td>rdintens</td>     <th>  R-squared:         </th> <td>   0.173</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.114</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   2.925</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th>  <td>0.0702</td> 
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -60.496</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>    31</td>      <th>  AIC:               </th> <td>   127.0</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>    28</td>      <th>  BIC:               </th> <td>   131.3</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable34">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    2.2969</td> <td>    0.592</td> <td>    3.881</td> <td> 0.001</td> <td>    1.085</td> <td>    3.509</td>
</tr>
<tr>
  <th>sales</th>     <td>    0.0002</td> <td> 8.42e-05</td> <td>    2.204</td> <td> 0.036</td> <td> 1.31e-05</td> <td>    0.000</td>
</tr>
<tr>
  <th>profmarg</th>  <td>    0.0478</td> <td>    0.044</td> <td>    1.075</td> <td> 0.291</td> <td>   -0.043</td> <td>    0.139</td>
</tr>
</table>
<table class="simpletable35">
<tr>
  <th>Omnibus:</th>       <td>19.377</td> <th>  Durbin-Watson:     </th> <td>   1.694</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  25.558</td>
</tr>
<tr>
  <th>Skew:</th>          <td> 1.623</td> <th>  Prob(JB):          </th> <td>2.82e-06</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 6.040</td> <th>  Cond. No.          </th> <td>8.56e+03</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 8.56e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems.

</div>

## Chapter 10: Basic Regression Analysis with Time Series Data

### **`Example 10.2:` Effects of Inflation and Deficits on Interest Rates**

$$\widehat{i3} = \beta_0 + \beta_1inf_t + \beta_2def_t$$
Data from the _Economic Report of the President, 2004_, Tables B-64, B-73, and B-79.


```python
dataWoo("intdef", description=True)
```

    name of dataset: intdef
    no of variables: 13
    no of observations: 56
    
    +----------+----------------------------------+
    | variable | label                            |
    +----------+----------------------------------+
    | year     | 1948 to 2003                     |
    | i3       | 3 month T-bill rate              |
    | inf      | CPI inflation rate               |
    | rec      | federal receipts, % GDP          |
    | out      | federal outlays, % GDP           |
    | def      | out - rec                        |
    | i3_1     | i3[_n-1]                         |
    | inf_1    | inf[_n-1]                        |
    | def_1    | def[_n-1]                        |
    | ci3      | i3 - i3_1                        |
    | cinf     | inf - inf_1                      |
    | cdef     | def - def_1                      |
    | y77      | =1 if year >= 1977; change in FY |
    +----------+----------------------------------+
    
    Economic Report of the President, 2004, Tables B-64, B-73, and B-79.



```python
intdef = dataWoo("intdef")
```


```python
intdef['infla'] = intdef['inf']
intdef['defla'] = intdef['def']
```


```python
tbill_model = smf.ols("i3 ~ infla + defla", data=intdef).fit()
tbill_model.summary()
```


<div>

<table class="simpletable36">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>           <td>i3</td>        <th>  R-squared:         </th> <td>   0.602</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.587</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   40.09</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>2.48e-11</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -112.16</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>    56</td>      <th>  AIC:               </th> <td>   230.3</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>    53</td>      <th>  BIC:               </th> <td>   236.4</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable37">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>    1.7333</td> <td>    0.432</td> <td>    4.012</td> <td> 0.000</td> <td>    0.867</td> <td>    2.600</td>
</tr>
<tr>
  <th>infla</th>     <td>    0.6059</td> <td>    0.082</td> <td>    7.376</td> <td> 0.000</td> <td>    0.441</td> <td>    0.771</td>
</tr>
<tr>
  <th>defla</th>     <td>    0.5131</td> <td>    0.118</td> <td>    4.334</td> <td> 0.000</td> <td>    0.276</td> <td>    0.751</td>
</tr>
</table>
<table class="simpletable38">
<tr>
  <th>Omnibus:</th>       <td> 0.260</td> <th>  Durbin-Watson:     </th> <td>   0.716</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.878</td> <th>  Jarque-Bera (JB):  </th> <td>   0.015</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.028</td> <th>  Prob(JB):          </th> <td>   0.992</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.058</td> <th>  Cond. No.          </th> <td>    9.28</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

</div>

### **`Example 10.11:` Seasonal Effects of Antidumping Filings**

C.M. Krupp and P.S. Pollard (1999), _Market Responses to Antidumpting Laws: Some Evidence from the U.S. Chemical Industry_, Canadian Journal of Economics 29, 199-227. Dr. Krupp kindly provided the data. They are monthly data covering February 1978 through December 1988. 


```python
dataWoo("barium", description=True)
```

    name of dataset: barium
    no of variables: 31
    no of observations: 131
    
    +----------+---------------------------------+
    | variable | label                           |
    +----------+---------------------------------+
    | chnimp   | Chinese imports, bar. chl.      |
    | bchlimp  | total imports bar. chl.         |
    | befile6  | =1 for all 6 mos before filing  |
    | affile6  | =1 for all 6 mos after filing   |
    | afdec6   | =1 for all 6 mos after decision |
    | befile12 | =1 all 12 mos before filing     |
    | affile12 | =1 all 12 mos after filing      |
    | afdec12  | =1 all 12 mos after decision    |
    | chempi   | chemical production index       |
    | gas      | gasoline production             |
    | rtwex    | exchange rate index             |
    | spr      | =1 for spring months            |
    | sum      | =1 for summer months            |
    | fall     | =1 for fall months              |
    | lchnimp  | log(chnimp)                     |
    | lgas     | log(gas)                        |
    | lrtwex   | log(rtwex)                      |
    | lchempi  | log(chempi)                     |
    | t        | time trend                      |
    | feb      | =1 if month is feb              |
    | mar      | =1 if month is march            |
    | apr      |                                 |
    | may      |                                 |
    | jun      |                                 |
    | jul      |                                 |
    | aug      |                                 |
    | sep      |                                 |
    | oct      |                                 |
    | nov      |                                 |
    | dec      |                                 |
    | percchn  | % imports from china            |
    +----------+---------------------------------+
    
    C.M. Krupp and P.S. Pollard (1999), "Market Responses to Antidumpting
    Laws: Some Evidence from the U.S. Chemical Industry," Canadian Journal
    of Economics 29, 199-227. Dr. Krupp kindly provided the data. They are
    monthly data covering February 1978 through December 1988.



```python
barium = dataWoo("barium")

barium_imports = smf.ols(
    "lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6",
    data=barium).fit()
barium_imports.summary()
```


<div>

<table class="simpletable39">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lchnimp</td>     <th>  R-squared:         </th> <td>   0.305</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.271</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   9.064</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>3.25e-08</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -114.79</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   131</td>      <th>  AIC:               </th> <td>   243.6</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   124</td>      <th>  BIC:               </th> <td>   263.7</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     6</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable40">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>  -17.8030</td> <td>   21.045</td> <td>   -0.846</td> <td> 0.399</td> <td>  -59.458</td> <td>   23.852</td>
</tr>
<tr>
  <th>lchempi</th>   <td>    3.1172</td> <td>    0.479</td> <td>    6.505</td> <td> 0.000</td> <td>    2.169</td> <td>    4.066</td>
</tr>
<tr>
  <th>lgas</th>      <td>    0.1964</td> <td>    0.907</td> <td>    0.217</td> <td> 0.829</td> <td>   -1.598</td> <td>    1.991</td>
</tr>
<tr>
  <th>lrtwex</th>    <td>    0.9830</td> <td>    0.400</td> <td>    2.457</td> <td> 0.015</td> <td>    0.191</td> <td>    1.775</td>
</tr>
<tr>
  <th>befile6</th>   <td>    0.0596</td> <td>    0.261</td> <td>    0.228</td> <td> 0.820</td> <td>   -0.457</td> <td>    0.576</td>
</tr>
<tr>
  <th>affile6</th>   <td>   -0.0324</td> <td>    0.264</td> <td>   -0.123</td> <td> 0.903</td> <td>   -0.556</td> <td>    0.491</td>
</tr>
<tr>
  <th>afdec6</th>    <td>   -0.5652</td> <td>    0.286</td> <td>   -1.978</td> <td> 0.050</td> <td>   -1.131</td> <td>    0.001</td>
</tr>
</table>
<table class="simpletable41">
<tr>
  <th>Omnibus:</th>       <td> 9.160</td> <th>  Durbin-Watson:     </th> <td>   1.458</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.010</td> <th>  Jarque-Bera (JB):  </th> <td>   9.978</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.491</td> <th>  Prob(JB):          </th> <td> 0.00681</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.930</td> <th>  Cond. No.          </th> <td>9.62e+03</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 9.62e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems.

</div>

Estimate a new model, `barium_seasonal` which accounts for seasonality by adding dummy variables contained in the data. 


```python
barium_seasonal = smf.ols(
    "lchnimp ~ lchempi + lgas + lrtwex + befile6 + affile6 + afdec6 + \
    feb + mar + apr + may + jun + jul + aug + sep + oct + nov + dec",
    data=barium).fit()
barium_seasonal.summary()
```


<div>

<table class="simpletable42">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>lchnimp</td>     <th>  R-squared:         </th> <td>   0.358</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.262</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   3.712</td>
</tr>
<tr>
  <th>Date:</th>             <td>Sat, 09 May 2020</td> <th>  Prob (F-statistic):</th> <td>1.28e-05</td>
</tr>
<tr>
  <th>Time:</th>                 <td>15:36:56</td>     <th>  Log-Likelihood:    </th> <td> -109.54</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   131</td>      <th>  AIC:               </th> <td>   255.1</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   113</td>      <th>  BIC:               </th> <td>   306.8</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>    17</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable43">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>   16.7788</td> <td>   32.429</td> <td>    0.517</td> <td> 0.606</td> <td>  -47.468</td> <td>   81.026</td>
</tr>
<tr>
  <th>lchempi</th>   <td>    3.2651</td> <td>    0.493</td> <td>    6.624</td> <td> 0.000</td> <td>    2.288</td> <td>    4.242</td>
</tr>
<tr>
  <th>lgas</th>      <td>   -1.2781</td> <td>    1.389</td> <td>   -0.920</td> <td> 0.359</td> <td>   -4.030</td> <td>    1.474</td>
</tr>
<tr>
  <th>lrtwex</th>    <td>    0.6630</td> <td>    0.471</td> <td>    1.407</td> <td> 0.162</td> <td>   -0.271</td> <td>    1.597</td>
</tr>
<tr>
  <th>befile6</th>   <td>    0.1397</td> <td>    0.267</td> <td>    0.524</td> <td> 0.602</td> <td>   -0.389</td> <td>    0.668</td>
</tr>
<tr>
  <th>affile6</th>   <td>    0.0126</td> <td>    0.279</td> <td>    0.045</td> <td> 0.964</td> <td>   -0.539</td> <td>    0.565</td>
</tr>
<tr>
  <th>afdec6</th>    <td>   -0.5213</td> <td>    0.302</td> <td>   -1.726</td> <td> 0.087</td> <td>   -1.120</td> <td>    0.077</td>
</tr>
<tr>
  <th>feb</th>       <td>   -0.4177</td> <td>    0.304</td> <td>   -1.372</td> <td> 0.173</td> <td>   -1.021</td> <td>    0.185</td>
</tr>
<tr>
  <th>mar</th>       <td>    0.0591</td> <td>    0.265</td> <td>    0.223</td> <td> 0.824</td> <td>   -0.465</td> <td>    0.584</td>
</tr>
<tr>
  <th>apr</th>       <td>   -0.4515</td> <td>    0.268</td> <td>   -1.682</td> <td> 0.095</td> <td>   -0.983</td> <td>    0.080</td>
</tr>
<tr>
  <th>may</th>       <td>    0.0333</td> <td>    0.269</td> <td>    0.124</td> <td> 0.902</td> <td>   -0.500</td> <td>    0.567</td>
</tr>
<tr>
  <th>jun</th>       <td>   -0.2063</td> <td>    0.269</td> <td>   -0.766</td> <td> 0.445</td> <td>   -0.740</td> <td>    0.327</td>
</tr>
<tr>
  <th>jul</th>       <td>    0.0038</td> <td>    0.279</td> <td>    0.014</td> <td> 0.989</td> <td>   -0.548</td> <td>    0.556</td>
</tr>
<tr>
  <th>aug</th>       <td>   -0.1571</td> <td>    0.278</td> <td>   -0.565</td> <td> 0.573</td> <td>   -0.708</td> <td>    0.394</td>
</tr>
<tr>
  <th>sep</th>       <td>   -0.1342</td> <td>    0.268</td> <td>   -0.501</td> <td> 0.617</td> <td>   -0.664</td> <td>    0.396</td>
</tr>
<tr>
  <th>oct</th>       <td>    0.0517</td> <td>    0.267</td> <td>    0.194</td> <td> 0.847</td> <td>   -0.477</td> <td>    0.580</td>
</tr>
<tr>
  <th>nov</th>       <td>   -0.2463</td> <td>    0.263</td> <td>   -0.937</td> <td> 0.351</td> <td>   -0.767</td> <td>    0.274</td>
</tr>
<tr>
  <th>dec</th>       <td>    0.1328</td> <td>    0.271</td> <td>    0.489</td> <td> 0.626</td> <td>   -0.405</td> <td>    0.671</td>
</tr>
</table>
<table class="simpletable44">
<tr>
  <th>Omnibus:</th>       <td> 9.169</td> <th>  Durbin-Watson:     </th> <td>   1.325</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.010</td> <th>  Jarque-Bera (JB):  </th> <td>   9.324</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.540</td> <th>  Prob(JB):          </th> <td> 0.00945</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.736</td> <th>  Cond. No.          </th> <td>1.47e+04</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.47e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems.

</div>

Now, compute the `anova` between the two models.


```python
from statsmodels.stats.api import anova_lm

anova_lm(barium_imports, barium_seasonal)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>df_resid</th>
      <th>ssr</th>
      <th>df_diff</th>
      <th>ss_diff</th>
      <th>F</th>
      <th>Pr(&gt;F)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>0</td>
      <td>124.0</td>
      <td>44.247088</td>
      <td>0.0</td>
      <td>NaN</td>
      <td>NaN</td>
      <td>NaN</td>
    </tr>
    <tr>
      <td>1</td>
      <td>113.0</td>
      <td>40.843896</td>
      <td>11.0</td>
      <td>3.403192</td>
      <td>0.855943</td>
      <td>0.585201</td>
    </tr>
  </tbody>
</table>
</div>


# おわりに
- Qiitaでも[github](https://github.com/lystahi/wooldridge_python_notebook/blob/master/Introductory-Econometrics-Examples.ipynb)でもご指摘やコメントなどしていいただけると嬉しいです。
- [ツイッター](https://twitter.com/lystahi) もやってます。