<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Statsmodels v.s. Scikit-learn
 Ruixue Gong, NYU
This notebook helps economists better understand the differences and similarities between machine learning package Scikit-learn and traditional statistical package Statsmodels. 
First of all, we provide an overview of Scikit-learn and Statsmodels. Secondly, we make a brief contrast of these two packages. At last, we talk about the applications on several topics such as linear regression models, logistic models and time series analysis with the use of Statsmodels and Scikit-learn respectively.
Overview of Scikit-learn
Step1:  Overview of Statsmodels
Linear regression models
Generalized linear models
Discrete choice models
Robust linear models
Many models and functions for time series analysis
Nonparametric estimators
A collection of datasets for examples
A wide range of statistical tests
Input-output tools for producing tables in a number of formats (Text, LaTex, HTML) and for reading Stata files into NumPy and Pandas.
Plotting functions
Extensive unit tests to ensure correctness of results
Contrast
Statsmodels
traditional model where we want to know how well a given model fits the data, and what variables "explain" or affect the outcome, or what the size of the effect is
analyzing the training data including hypothesis tests and goodness-of-fit measures
Scikit-learn
machine learning algorithms where the main supported task is choosing the "best" model for prediction
model selection for out-of-sample prediction and hence cross-validation on "test data".
Contents
Generalized Linear Models
OLS
Testing
Logistic Regression
Time Series Analysis
Step2:  OLS and Testing
Linear regression is in its basic form the same in statsmodels and in scikit-learn. However, Statsmodels can generate a summary of OLS result like Stata and R. Lecture 11 had already given an introduction to regression analysis in Sciki-learn, so we only present the different features owned by Statsmodels here.
Now consider the real model as follows
Step3:  We can get the summary of the regression result just like R.
Summary
results.summary()
We can also extract the values that we are interested in.
Coefficients
Step4:  Joint Test
We can test $$R\beta = 0$$ or formula-like $$\beta_2 = \beta_3 = 0$$
Step5:  Maximum Likelihood Estimation
Statsmodel has many built-in regression models using maximum likelihood estimation. However, Scikit-learn doesn't.
Step6:  Now I want to compare binary predictions.
Logistic Regression using Statsmodels.
Supporting Vector Machine (SVM) with linear kernal using Scikit-learn
I use the spector dataset which contained in Statsmodels package. To visualize the results, I set two exogenous variables, x1 and x2.
We all have familiarized with the mathematical formation of logistic regression. So I put a brief introduction graph about supporting vector machine algorithm as follows.
Step7:  Time Series Analysis
As I mentioned in the beginning, Scikit-learn focus on prediction of time series by machine learning algorithm. It doesn't contain traditional time series analysis such as ARMA, VAR, and relevant statictical tests. But you can find these in Statsmodels package.
<ASSISTANT_TASK:>
Python Code:
from IPython.display import Image
Image(filename='scikit-learn-flow-chart.jpg')  #source: web
Explanation: Statsmodels v.s. Scikit-learn
 Ruixue Gong, NYU
This notebook helps economists better understand the differences and similarities between machine learning package Scikit-learn and traditional statistical package Statsmodels. 
First of all, we provide an overview of Scikit-learn and Statsmodels. Secondly, we make a brief contrast of these two packages. At last, we talk about the applications on several topics such as linear regression models, logistic models and time series analysis with the use of Statsmodels and Scikit-learn respectively.
Overview of Scikit-learn
End of explanation
%matplotlib inline
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.sandbox.regression.predstd import wls_prediction_std
from scipy import stats
import pandas as pd
from sklearn import svm
from sklearn.linear_model import LogisticRegression
from statsmodels.base.model import GenericLikelihoodModel
Explanation: Overview of Statsmodels
Linear regression models
Generalized linear models
Discrete choice models
Robust linear models
Many models and functions for time series analysis
Nonparametric estimators
A collection of datasets for examples
A wide range of statistical tests
Input-output tools for producing tables in a number of formats (Text, LaTex, HTML) and for reading Stata files into NumPy and Pandas.
Plotting functions
Extensive unit tests to ensure correctness of results
Contrast
Statsmodels
traditional model where we want to know how well a given model fits the data, and what variables "explain" or affect the outcome, or what the size of the effect is
analyzing the training data including hypothesis tests and goodness-of-fit measures
Scikit-learn
machine learning algorithms where the main supported task is choosing the "best" model for prediction
model selection for out-of-sample prediction and hence cross-validation on "test data".
Contents
Generalized Linear Models
OLS
Testing
Logistic Regression
Time Series Analysis
End of explanation
n = 20
x = np.linspace(0, 5, n)
sigma = 0.3
beta = np.array([1, 0.5, -0.02,5])    # real coefficient
e = np.random.normal(size=n)
X = np.column_stack((x, np.sin(x), (x-3)**2, np.ones(n))) 
y_true = np.dot(X,beta)
y = y_true + e
#do regression
model = sm.OLS(y, X)   #Pick a class. GLS, WLS...
results = model.fit()
Explanation: OLS and Testing
Linear regression is in its basic form the same in statsmodels and in scikit-learn. However, Statsmodels can generate a summary of OLS result like Stata and R. Lecture 11 had already given an introduction to regression analysis in Sciki-learn, so we only present the different features owned by Statsmodels here.
Now consider the real model as follows: 
$$y = \beta_1 x + \beta_2 sin(x) + \beta_3 (x-3)^2 + \beta_0, \quad \text{ where } \epsilon \sim N(0,1)$$
End of explanation
print(results.summary())
print('Coefficients: ', results.params)
print('R2: ', results.rsquared)
print('Standard errors: ', results.bse)
print('Fitted values: ', results.fittedvalues)
print('Predicted values: ', results.predict())
prstd, iv_l, iv_u = wls_prediction_std(results)
#wls_prediction_std returns standard deviation and 
#confidence interval of my fitted model data
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, y_true, 'b-', label="Real")
ax.plot(x, results.fittedvalues, 'r--.', label="OLS")
ax.plot(x, iv_u, 'r--')
ax.plot(x, iv_l, 'r--')
ax.legend(loc='best')
Explanation: We can get the summary of the regression result just like R.
Summary
results.summary()
We can also extract the values that we are interested in.
Coefficients:        results.params
R square: results.rsquared
Fitted values: results.fittedvalues
Predicted values:   results.predict()
Standard errors of each coefficient:  results.bse
End of explanation
#test beta_2 = beta_3 = 0
print(results.f_test("x2 = x3 = 0"))
#test R beta = 0
R = [[0, 1, 0, 0], [0, 0, 1, 0]]
print(results.f_test(R))
Explanation: Joint Test
We can test $$R\beta = 0$$ or formula-like $$\beta_2 = \beta_3 = 0$$
End of explanation
data = sm.datasets.spector.load_pandas()  #use the dataset spector
exog = data.exog
endog = data.endog
print(sm.datasets.spector.NOTE)
print(data.exog.head())
exog1 = sm.add_constant(exog, prepend=True)   #combine X matrix with constant
#plug in the log-likelihood function of my own model
class MyProbit(GenericLikelihoodModel):
    def loglike(self, params):
        exog = self.exog
        endog = self.endog
        q = 2 * endog - 1
        y = stats.norm.logcdf(q*np.dot(exog, params)).sum()   
        return y
sm_probit_manual = MyProbit(endog, exog1).fit()
print(sm_probit_manual.summary())
Explanation: Maximum Likelihood Estimation
Statsmodel has many built-in regression models using maximum likelihood estimation. However, Scikit-learn doesn't.
End of explanation
Image(filename='svm.jpg' )     #source : web
X = np.delete(data.exog.as_matrix(),2,1)
y = np.asarray(data.endog.tolist())
#classifiers
clf = svm.SVC(kernel='linear', C = 1)  
#fitting by Scikit-learn
model_svm = clf.fit(X,y)
#fitting by Statsmodels
XX = np.insert(X,0,1,axis = 1) #add a column with contant 1
model_logit_stats = sm.Logit(y,XX)
res = model_logit_stats.fit()
w = clf.coef_[0]
a = -w[0] / w[1]
c = -res.params[1] / res.params[2]
x1 = np.linspace(2,4)
x2 = a * x1 - clf.intercept_[0] / w[1]     #svm 
x3 = c * x1 - res.params[0] / res.params[2]     #logistic
h1 = plt.plot(x1,x2,'k-', color = 'green', label = 'SVM')
h2 = plt.plot(x1,x3,'k-', color = 'blue', label = 'Logit Statsmodels')
plt.scatter(X[:, 0], X[:,1], c = y, alpha = 0.4)
plt.legend()
plt.title('Comparision: LogisticRegression and SVM')
plt.show()
Explanation: Now I want to compare binary predictions.
Logistic Regression using Statsmodels.
Supporting Vector Machine (SVM) with linear kernal using Scikit-learn
I use the spector dataset which contained in Statsmodels package. To visualize the results, I set two exogenous variables, x1 and x2.
We all have familiarized with the mathematical formation of logistic regression. So I put a brief introduction graph about supporting vector machine algorithm as follows.
End of explanation
dta = sm.datasets.sunspots.load_pandas().data  #use the dataset sunspot
print(sm.datasets.sunspots.NOTE)
dta.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008')) #set time series
del dta["YEAR"]
dta.plot(figsize=(12,8))
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(dta.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(dta, lags=40, ax=ax2)
Explanation: Time Series Analysis
As I mentioned in the beginning, Scikit-learn focus on prediction of time series by machine learning algorithm. It doesn't contain traditional time series analysis such as ARMA, VAR, and relevant statictical tests. But you can find these in Statsmodels package.
End of explanation
<END_TASK>