<jupyter_start><jupyter_text>Introduction into Data Science with Python  PrefaceOver the last decades statistical analysis becomes more and more complex. Statisticians needs for software and machines to compute large matrices with complex operations arises. Since measurement has become more accurate, storage is easy to access/ share and the awareness that tracked data can be useful later for researcher of all fields, the amount of data available is tremendous. Therefore, modern statistical software also needs handle with large dataset, while providing efficient algorithms.Statistics Software available:    A. Proprietary        - SPSS        - STATA        - GRETL        - EViews        - GAUSS        - Matlab        ...    B: Open Source        - R        - Python        - Julia        - Octave        ...        Currently, *Open Source* software is very popular:    - free to use    - support in forums from statisticians, scientists and professionals     - infinitely expandable by user-based packages    - very easy to start group project since it's free    Why to use **Python**?R is widely used in all fields, providing a large amount on packages and with a very responsive fan-base providing support. It's still very slow in computing, requires a lot of computing performance and is not as easy to read the code as Python. For example, R-code requires many different *brackets* depending on their usage. Mixing up these *brackets* or missing one leads to errors in the code.Here Python shines with code easy to read and as a programming language mainly used to engineer software it provides a huge potential. Here the first difference appears. Since, Python is a programming language we first need to load packages which include functions and operations we later want to use for our analysis. For users of proprietary statistic software it might be uncontentious, to load packages to compute basic statistics.  We use Python together with the Jupyter-Notebook Editor. Further information regarding the Editor can be found [here](https://www.dataquest.io/blog/jupyter-notebook-tips-tricks-shortcuts/)<jupyter_code>###############################################
### Introduction into Python
###############################################

## Preface - loading packages
import sys                        # package for System-specific parameters and functions in Python
import numpy as np                # scientific operations
import pandas as pd               # package for DataFrames and related functions
import statsmodels.api as sm      # selecting *.api since there are different versions available (*.api is designed to work with pandas)
import statsmodels.formula.api as smf  # for more advanced regression models
from patsy import dmatrices       # using the dmatrices function from the patsy package
import matplotlib.pyplot as plt   # package to producing graphs and plots
import scipy.stats as ss          # scientific operations
import seaborn as sns             # package based on matplotlib

## Preface - setting some global options
# %matplotlib inline              # procudes the output inline with the Notebook<jupyter_output><empty_output><jupyter_text>Basic operations in Python<jupyter_code>## Basic operations
# storing and computation
a=4
b=5
a+b
# Output 
print("Hello")
# IF Statements
if a == 4:
    print('a is four')
elif b == 5:
    print('b is five')
else:
    print('a is not four and b not five')
# Loops
for i in ["First", "Second", "Thrid"]:
    print(i, "rank")
for i in range(1,4): 
    print(i)
my_list = [1, 2, 3]
for i in my_list: 
    print(i)
# Functions
my_list = [0,1,2,3,4,5,6,7,8,9,10]
def func(x):
    res=0
    for i in x:
        if i>5:
         res = res + i
    return res

func(my_list)
# consulting the help-file
help(sum)<jupyter_output>Help on built-in function sum in module builtins:

sum(iterable, start=0, /)
    Return the sum of a 'start' value (default: 0) plus an iterable of numbers
    
    When the iterable is empty, return the start value.
    This function is intended specifically for use with numeric values and may
    reject non-numeric types.<jupyter_text>Importing Data using Pandas package<jupyter_code># reading a stata file using the pandas package as "pd"
df = pd.read_stata('nswre74.dta')
df.head()
df.describe()
# producing a simple scatter plot
plt.scatter(df.re75, df.re74)
plt.xlabel('re75')
plt.ylabel('re74')<jupyter_output><empty_output><jupyter_text>Regressions in Python using Statsmodels package<jupyter_code># using statsmodels as "sm"
y, x = dmatrices('re78 ~ treat + black', data = df)

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

print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                   re78   R-squared:                       0.032
Model:                            OLS   Adj. R-squared:                  0.027
Method:                 Least Squares   F-statistic:                     7.239
Date:                Wed, 17 Apr 2019   Prob (F-statistic):           0.000806
Time:                        13:57:46   Log-Likelihood:                -4539.6
No. Observations:                 445   AIC:                             9085.
Df Residuals:                     442   BIC:                             9097.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Robust OLS<jupyter_code># Robust Regression with FE
y, x = dmatrices('re78 ~ treat + black', data = df)

model_spec = sm.RLM(y, x, M=sm.robust.norms.HuberT())
results = model_spec.fit()

print(results.summary())<jupyter_output>Robust linear Model Regression Results                    
==============================================================================
Dep. Variable:                   re78   No. Observations:                  445
Model:                            RLM   Df Residuals:                      442
Method:                          IRLS   Df Model:                            2
Norm:                          HuberT                                         
Scale Est.:                       mad                                         
Cov Type:                          H1                                         
Date:                Wed, 17 Apr 2019                                         
Time:                        13:57:47                                         
No. Iterations:                     8                                         
==============================================================================
                 coef    std err          z      P>|[...]<jupyter_text>OLS with HAC standard errors<jupyter_code># Robust Regression with FE
y, x = dmatrices('re78 ~ treat + black', data = df)

model_spec = sm.OLS(y, x)
results = model_spec.fit(cov_type='HAC',cov_kwds={'maxlags':1})

print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                   re78   R-squared:                       0.032
Model:                            OLS   Adj. R-squared:                  0.027
Method:                 Least Squares   F-statistic:                     8.964
Date:                Wed, 17 Apr 2019   Prob (F-statistic):           0.000153
Time:                        13:57:48   Log-Likelihood:                -4539.6
No. Observations:                 445   AIC:                             9085.
Df Residuals:                     442   BIC:                             9097.
Df Model:                           2                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|[...]<jupyter_text>Random numbers and data generating processes (DGP)  Random numbers<jupyter_code>np.random.seed(123)       # set seed to ensure reproducibility
A = np.random.rand(3,2)   # generating a 3X2 random matrix 
A
B = np.random.binomial(100, 0.5)   # generating binomial(n,p)
B
C = np.random.uniform(0, 1, 100)   # generating binomial(low,upper, n)
C
# Creating random data, calculating density function and plot
figure, aux = plt.subplots(1, 1)
x = np.linspace(-5, 5, 1500)      # Return evenly spaced numbers over a specified interval. (start, stop, number)
pdf = ss.norm.pdf(x, 0, 1)       # normal distributed probability function with zero mean and variance of 1
aux.plot(x, pdf)

# plotting options
aux.set_xlabel('realizations')
aux.set_ylabel('Density')
#aux.set_xticklabels(['', '', '', 0.5, '', '', ''])
#aux.set_xlim([-3, 3])
#aux.set_ylim([0, 0.5]);<jupyter_output><empty_output><jupyter_text>Data generating process (DGP)<jupyter_code>Num= 1000       # number of observations

alpha = 0.8
beta  = 0.5

data = np.tile(np.nan, (Num, 2))
for i in range(Num):
    X      = np.random.uniform() 
    res    = np.random.normal()
    Y      = alpha + beta*X + res
    data[i, :] = [Y, X]

df=pd.DataFrame(data, columns=['Y','X'])

y, x = dmatrices('Y ~ X', data = df)

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

print(results.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.012
Model:                            OLS   Adj. R-squared:                  0.011
Method:                 Least Squares   F-statistic:                     12.55
Date:                Wed, 17 Apr 2019   Prob (F-statistic):           0.000415
Time:                        14:52:39   Log-Likelihood:                -1427.3
No. Observations:                1000   AIC:                             2859.
Df Residuals:                     998   BIC:                             2868.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Resampling data<jupyter_code>df.describe()    #Lalonde Dataset
# creating variables for further use
is_treated = df['treat'] == 1
df['Y'] = df['re78']
df['Y_0'] = df.loc[~is_treated, 're78']
df['Y_1'] = df.loc[is_treated, 're78']
df['D'] = np.nan
df.loc[~is_treated, 'D'] = 0
df.loc[is_treated, 'D'] = 1

# set the seed for reproduceability
np.random.seed(123)
# caluclate the size of Treatment group
status_counts = df['D'].value_counts().to_dict()

stats = list()
for _ in range(100):
    Y_1_sampled = df['Y_1'].dropna().sample(n=status_counts[0], replace=True).values
    Y_0_sampled = df['Y_0'].dropna().sample(n=status_counts[1], replace=True).values

    df_boot = df.copy()

    df_boot.loc[df_boot['D'] == 1, 'Y_0'] = Y_0_sampled
    df_boot.loc[df_boot['D'] == 0, 'Y_1'] = Y_1_sampled
    
    stat = (df_boot['Y_1'] - df_boot['Y_0']).mean()
    stats.append(stat)
print('Our estimated treatment effect is {:10.2f}'.format(np.mean(stats)))
fig, ax = plt.subplots(1, 1)
ax.hist(stats)
ax.set_ylabel('Statistic')
ax.set_xlabel('Frequency');<jupyter_output>Our estimated treatment effect is    1787.05<jupyter_text>Robust Regression example<jupyter_code>np.random.seed(123) 
nsample = 50
x1 = np.linspace(0, 20, nsample)
X = np.column_stack((x1, (x1-5)**2))
X = sm.add_constant(X)
sig = 0.3                                         # smaller error variance makes OLS<->RLM contrast bigger
beta = [5, 0.5, -0.0]
y_true2 = np.dot(X, beta)
y2 = y_true2 + sig*1. * np.random.normal(size=nsample)
y2[[39,41,43,45,48]] -= 5                         # add some outliers (10% of nsample)

res = sm.OLS(y2, X).fit()
resrlm = sm.RLM(y2, X).fit()

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(x1, y2, 'o',label="data")
ax.plot(x1, y_true2, 'b-', label="True")
ax.plot(x1, res.fittedvalues, 'r-', label="OLS")
ax.plot(x1, resrlm.fittedvalues, 'g.-', label="RLM")
ax.legend(loc="best")<jupyter_output><empty_output><jupyter_text>Loading Data using API  World bank data<jupyter_code>#!pip install wbdata # we need to install a third-party package
import wbdata        # easy access to World Bank data
 
# source: https://data.worldbank.org/indicator/sp.dyn.cbrt.in
indicators = {'SP.DYN.CBRT.IN':'Birth rate'}
df = wbdata.get_dataframe(indicators, country=["JP"], convert_date=False)

df.sort_index(inplace=True)       # sort by year
df.index = df.index.astype(int)   # format index to integer type
df.head()<jupyter_output><empty_output><jupyter_text>Quandl API for multiple sources<jupyter_code>#!pip install quandl
# note for more than 50 data requests you need to sign up
# quandl.ApiConfig.api_key = "YOUR_KEY_HERE"

import quandl

df = quandl.get("FRED/GDP", collapse="monthly")
df.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14729209154769998, 'fraction_numerical': 0.027418989349648766, 'mean_word_length': 2.240332843857073, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 49, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '178568', 'license': 'MIT', 'n_tokens_mistral': 4054, 'n_tokens_neox': 3704, 'n_words': 1159, 'path': 'Introduction_into_Python.ipynb', 'repo_name': 'fbalensiefer/Introduction_Python', 'size': '13239'}