#!/usr/bin/env python
# coding: utf-8

# In[28]:


import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy as sp
import scipy.stats
import statsmodels.api as sm


# # <font face="gotham" color="purple"> ANOVA </font>

# If you have studied statistics, you certainly know the famous **Analysis of Variance** (ANOVA), you can skip this section, but if you haven't, read on.
# 
# Simply speaking, the ANOVA is a technique of comparing means of multiple$(\geq 3)$ populations, the name derives from the way how calculations are performed. 
# 
# For example, a common hypotheses of ANOVA are 
# $$
# H_0:\quad \mu_1=\mu_2=\mu_3=\cdots=\mu_n\\
# H_1:\quad \text{At least two means differ}
# $$
# In order to construct $F$-statistic, we need to introduce two more statistics, the first one is **Mean Square for Treatments** (MST), $\bar{\bar{x}}$ is the grand mean, $\bar{x}_i$ is the sample mean of sample $x_i$, $n_i$ is the number of observable in sample $i$
# $$
# MST=\frac{SST}{k-1},\qquad\text{where } SST=\sum_{i=1}^kn_i(\bar{x}_i-\bar{\bar{x}})^2
# $$
# And the second one is **Mean Square for Error** (MSE), $s_i$ is the sample variance of sample $i$
# $$
# MSE=\frac{SSE}{n-k},\qquad\text{where } SSE =(n_1-1)s_1^2+(n_2-1)s_2^2+\cdots+(n_k-1)s_k^2
# $$
# Join them together, an $F$-statistic is constructed
# $$
# F=\frac{MST}{MSE}
# $$
# If the $F$-statistic is larger than critical value with its corresponding degree of freedom, we reject null hypothesis.

# # <font face="gotham" color="purple"> Dummy Variable </font>

# Here's dataset with dummy variables, which are either $1$ or $0$. 

# In[29]:


df = pd.read_excel("Basic_Econometrics_practice_data.xlsx", sheet_name="Hight_ANOVA")


# In[30]:


df


# The dataset has five columns, the first column $Height$ is a sample of $88$ male height, other columns are dummy variables indication its qualitative feature, here is the nationality.
# 
# There are $4$ countries in the sample, Japan, Netherlands, Denmark and Finland, however there are only $3$ dummies in the data set, this is to avoid _perfect multicollinearity_, which is also called **dummy variable trap**, because the height data is the perfect linear combination of four dummy variables. 
# 
# If we use the model with only dummy variables as independent variable, we basically regressing a ANOVA model, i.e.
# $$
# Y_{i}=\beta_{1}+\beta_{2} D_{2 i}+\beta_{3 i} \mathrm{D}_{3 i}+\beta_{3 i} \mathrm{D}_{3 i}+u_{i}
# $$

# where $Y_i =$ the male height, $D_{2i}=1$ if the male is from Netherlands,  $D_{3i}=1$ if the male is from Denmark and $D_{4i}=1$ if the male is from Finland. Japan doesn't have a dummy variable, so we are using it as reference, which will be clearer later.
# 
# Now we run the regression and print the result. And how do we interpret the estimated coefficients?

# In[31]:


X = df[["NL_dummpy", "DM_dummpy", "FI_dummy"]]
Y = df["Height"]

X = sm.add_constant(X)  # adding a constant

model = sm.OLS(Y, X).fit()
print_model = model.summary()
print(print_model)


# First, all the $p$-value are significantly small, so our estimation is valid. Then we examine the coefficients one by one. 
# 
# The estimated constant $b_1 = 163.03$ is the mean height of Japanese male. The mean of Dutch male height is $b_1+b_2 = 163.03+17.71=180.74$, the mean of Danish male height is $b_1+b_3=163.03+12.21=175.24$, the mean of Finnish male height is $b_1+b_4=163.03+12.85=175.88$. 
# 
# As you can see, the Japanese male height is used as a _reference_, also called _base category_, rest are added onpon it. 

# This regression has the same effect as ANOVA test, all dummy coefficients are significant and so is $F$-statistic, which means we reject null height null hypothesis that all countries' male height are the same.

# # <font face="gotham" color="purple"> The ANCOVA Models </font>

# The example in the last section has only dummies in the independent variables, which is rare in practice. The more common situation is that independent variables have both quantitative and qualitative/dummy variables, and this is what we will do in this section. 
# 
# The model with both quantitative and qualitative variables are called **analysis of covariance** (ANCOVA) models. We have another dataset that contains the individual's parents' height. 

# In[54]:


df = pd.read_excel("Basic_Econometrics_practice_data.xlsx", sheet_name="Hight_ANCOVA")

X = df[["Height of Father", "Height of Mother", "NL_dummy", "DM_dummy", "FI_dummy"]]
Y = df["Height"]

X = sm.add_constant(X)  # adding a constant

model = sm.OLS(Y, X).fit()
print_model = model.summary()
print(print_model)


# In order to interpret the results, let's type the estimated model here
# $$
# \hat{Y}= 27.87+.33 X_{f} + .5 X_{m} + 5.36 D_{NL} + 2.90 D_{DM} + 1.02 D_{FI}
# $$
# where $X_{f}$ and $X_{m}$ are father's and mother's heights, $D$'s are dummy variables representing each country.

# This is actually a function to predict a male's height if you input parents height, for instance if we set $D_{NL} = D_{DM}= D_{FI}=0 $, the function of height of Japanese male is
# $$
# \hat{Y}= 27.87+.33 X_{f} + .5 X_{m}
# $$
# Or the function of Dutch male height with $D_{NL} = 1$ and $ D_{DM}= D_{FI}=0$
# $$
# \hat{Y}= 27.87+.33 X_{f} + .5 X_{m} + 5.36
# $$
# With these results, we can define Python functions to predict male height

# In[50]:


def jp_height(fh, mh):
    return 27.87 + fh * 0.33 + mh * 0.5


def nl_height(fh, mh):
    return 27.87 + fh * 0.33 + mh * 0.5 + 5.36


# A function to predict a Japanese male's height

# In[34]:


jp_height(175, 170)


# And function to predict a Dutch male's height

# In[35]:


nl_height(185, 185)


# # <font face="gotham" color="purple"> Slope Dummy Variables </font>

# What we have discussed above are all **intercept dummy variables** which means the dummy variable only change the intercept term, however, there are dummies which imposed on slope variables too. 
# 
# Back to the height example, what if we suspect that parents' height in NL could have more marginal effect on their sons' height, i.e. the model is
# $$
# Y= \beta_1 + \beta_2D_{NL} + (\beta_3+ \beta_4D_{NL})X_{f} +  (\beta_5+ \beta_6D_{NL})X_{m}+u
# $$
# Rewrite for estimation purpose
# $$
# Y= \beta_1 + \beta_2D_{NL} + \beta_3 X_f + \beta_4 D_{NL}X_f + \beta_5X_m + \beta_6 D_{NL}X_m+u
# $$

# Take a look at our data, we need to reconstruct it

# In[51]:


df.head()


# Drop the dummies of Denmark and Finland

# In[55]:


df_NL = df.drop(["DM_dummy", "FI_dummy"], axis=1)
df_NL.head()


# Also create the column of multiplication of $D_{NL} \cdot X_f$ and $D_{NL}\cdot X_m$, then construct independent variable matrix and dependent variable

# In[56]:


df_NL["D_NL_Xf"] = df_NL["Height of Father"] * df_NL["NL_dummy"]
df_NL["D_NL_Xm"] = df_NL["Height of Mother"] * df_NL["NL_dummy"]
X = df_NL[["NL_dummy", "Height of Father", "D_NL_Xf", "Height of Mother", "D_NL_Xm"]]
Y = df["Height"]


# In[59]:


X = sm.add_constant(X)  # adding a constant

model = sm.OLS(Y, X).fit()
print_model = model.summary()
print(print_model)


# Here's the estimated regression model
# $$
# \hat{Y}=  11.7747  + 120.9563D_{NL} + 0.3457 X_f  - 0.0946 D_{NL}X_f + 0.5903X_m  -0.5746 D_{NL}X_m
# $$
# If $D_{NL}=1$ then
# $$
# \hat{Y}=  132.731+0.2511X_f   -0.01X_m
# $$

# Again, we define a Python function to predict Dutch male height based on their parents' height

# In[69]:


def nl_height_marginal(fh, mh):
    return 132.731 + fh * 0.2511 + mh * 0.0157


# In[73]:


nl_height_marginal(185, 180)


# Prediction seem quite logical. 
# 
# However, as you can see from the results, the hypotheses test rejects our theory that Dutch parents could influence their sons' marginal height, i.e. coefficients of $D_{NL} \cdot X_f$ and $D_{NL}\cdot X_m$ fail to reject null hypothesis with $5\%$ significance level.

# In[ ]: