def inbound_forcast(target, exchange, geo, exchange_test, geo_test, submit, i):
    for col in col_list:
        # 宿泊者数のカラム名を指定
        target_col = col + suff

        target.index = range(0, 365)

        X = sm.add_constant(exchange, prepend=False)
        X_test = sm.add_constant(exchange_test, prepend=False)

        X.index = range(0,365)

        for g in range(0, len(target)):
            if target[target_col][g] == 0:
                target[target_col][g] = 1

        y = target[target_col].apply(np.log)

        model = sm.OLS(y, X)
        results = model.fit()
        print(results.summary())
        pred = results.predict()
        Y = y - pred
        L = len(Y)

        fftY = fft.fft(Y)
        freqs = fft.fftfreq(L)
        power = np.abs(fftY)
        phase = [np.arctan2(float(c.imag), float(c.real)) for c in fftY]

        wave = newwave_i(L, results, pred, power, freqs, phase, X_test)
        submit[i] = wave
        i += 1

    return submit, i

def error_rate_for_model(test_model, train_set, test_set, infer=False,
    infer_set=None):
    """Report error rate on test_doc sentiments, using
    supplied model and train_docs"""

    train_targets, train_regressors = \
        zip(*[( doc.sentiment, test_model.docvecs[doc.tags[0]] )
                for doc in train_set ])
    train_regressors = sm.add_constant(train_regressors)
    predictor = logistic_predictor_from_data(train_targets, train_regressors)

    test_data = test_set
    if infer:
        #if infer_subsample < 1.0:
        #    test_data = sample(test_data,
        #                       int(infer_subsample * len(test_data)))
        #test_regressors = [test_model.infer_vector(doc.words,
        #                       steps=infer_steps, alpha=infer_alpha)
        #                   for doc in test_data]
        test_data = [SentimentDocument(None, None, None, s)
                     for (v, s) in infer_set]
        test_regressors = [v for (v, s) in infer_set]
    else:
        test_regressors = [test_model.docvecs[doc.tags[0]]
                           for doc in test_set]
    test_regressors = sm.add_constant(test_regressors)

    # predict & evaluate
    test_predictions = predictor.predict(test_regressors)
    corrects = ( sum(np.rint(test_predictions) ==
                 [doc.sentiment for doc in test_data] ))
    errors = len(test_predictions) - corrects
    error_rate = float(errors) / len(test_predictions)
    return (error_rate, errors, len(test_predictions), predictor)

 def __init__(self):
   self.researched = util.source.read('F-F_Research_Data_5_Factors_2x3')
   self.portfolios = util.source.read('25_Portfolios_5x5')
   self.simpleFactor = sm.add_constant(self.researched.Mkt_RF)
   self.threeFactor = sm.add_constant(self.researched[['Mkt_RF', 'SMB', 'HML']])
   self.fourFactor = sm.add_constant(self.researched[['Mkt_RF', 'SMB', 'RMW', 'CMA']])
   self.fiveFactor = sm.add_constant(self.researched[['Mkt_RF', 'SMB', 'HML', 'RMW', 'CMA']])

def linear_model_plot(x_variable, y_variable):
    '''Function develops linear model for x and y variable inputs and plots regression line on top of scatter plot'''
    
    assert len(x_variable) > 1, 'length of x_variable should be larger than 1'
    assert len(y_variable) > 1, 'length of y_variable should be larger than 1'
    
    # assigning function variables to response and predictor variables
    y = y_variable # response variable
    X = x_variable # predictor variable
    X = sm.add_constant(X)  # Adds a constant term to the predictor (essential to obtain the constant in the formula)
    
    # Calculating the linear model for the two variables
    lm = sm.formula.OLS(y, X).fit()
    
    # Developing the plot of the linear model
    # making a range of the x variable to pass to the y prediction
    x_pred = np.linspace(x_variable.min(), x_variable.max())
    
    # Adding a constant to this range of x values (essential to obtain the constant in the formula)
    x_pred2 = sm.add_constant(x_pred)

    # Passing the linear model predictor the range of x values to model over
    y_pred = lm.predict(x_pred2)

    # Plotting these predicitons on the graph
    plt.plot(x_pred, y_pred, color='k', linewidth=2)

    # Obtaining linear regression 
    return plt.plot()

    def predict(self, test_X):
        dataset = self.__dataset
        intercept = self.__intercept
        XX_inv  = self.__XX_inv
        beta    = self.__beta
        
        train_X = sm.add_constant(dataset[:, :-1]) if intercept else dataset[:, :-1]
        test_X = sm.add_constant(vec(test_X)) if intercept else vec(test_X)
        train_Y = dataset[:, -1:]
        train_pred = np.dot(train_X, beta)
        
        # Confidence interval
        sig = (np.linalg.norm(train_Y-train_pred)**2/(train_X.shape[0]-train_X.shape[1]+1))**0.5
        s = []
        for row in range(test_X.shape[0]):
            x = test_X[[row], :]
            s.append(sig*(1 + np.dot(np.dot(x, XX_inv), x.T))**0.5)
            
        s = np.reshape(np.asarray(s), (test_X.shape[0], 1))

        test_pred = np.dot(test_X, beta)
        hi_ci = test_pred + 2*s
        lo_ci = test_pred - 2*s

        return test_pred, hi_ci, lo_ci

def overfit_stocks():
  # Load one year's worth of pricing data for five different assets
  start = datetime.date(1,1,2013)
  end = datetime.datetime(1,1,2014)
  x1 = get_pricing('PEP', )
  x2 = get_pricing('MCD', fields='price', start_date=start, end_date=end)
  x3 = get_pricing('ATHN', fields='price', start_date=start, end_date=end)
  x4 = get_pricing('DOW', fields='price', start_date=start, end_date=end)
  y = get_pricing('PG', fields='price', start_date=start, end_date=end)
  #
  # Build a linear model using only x1 to explain y
  slr = regression.linear_model.OLS(y, sm.add_constant(x1)).fit()
  slr_prediction = slr.params[0] + slr.params[1]*x1
  #
  # Run multiple linear regression using x1, x2, x3, x4 to explain y
  mlr = regression.linear_model.OLS(y, sm.add_constant(np.column_stack((x1,x2,x3,x4)))).fit()
  mlr_prediction = mlr.params[0] + mlr.params[1]*x1 + mlr.params[2]*x2 + mlr.params[3]*x3 + mlr.params[4]*x4
  #
  # Compute adjusted R-squared for the two different models
  print('SLR R-squared: %.5f' %slr.rsquared_adj)
  print('SLR p-value: %.5f' %slr.f_pvalue)
  print('MLR R-squared: %.5f' %mlr.rsquared_adj)
  print('MLR p-value: %.5f'  %mlr.f_pvalue)
  #
  # Plot y along with the two different predictions
  y.plot()
  slr_prediction.plot()
  mlr_prediction.plot()
  plt.ylabel('Price')
  plt.xlabel('Date')
  plt.legend(['PG', 'SLR', 'MLR']);

def GetCoef(start_train, end_train, StockReturns, CarhartDaily, SP500Returns, DataFolder):
    if os.path.isfile(r'%s\Coef_%s_%s.csv' % (DataFolder, start_train.date(), end_train.date())):
        Coef = pd.read_csv(r'%s\Coef_%s_%s.csv' % (DataFolder, start_train.date(), end_train.date()))
        return Coef
    else:
        Coef = pd.DataFrame()
        for ticker in StockReturns.ticker.unique():
            print "Getting regression coefficient for %s" % ticker
            tmpReturn = StockReturns[(StockReturns.ticker == ticker)]
            if not tmpReturn.empty:
                tmpData = tmpReturn.merge(CarhartDaily, left_on = 'endDate', right_on = 'date')
                tmpData = tmpData.merge(SP500Returns, on = 'endDate')
                tmpData['SP500-RF'] = tmpData['SP500Return']*100 - tmpData['RF']
                y = tmpData['return']*100 - tmpData['RF']
                X1 = tmpData[['Mkt-RF', 'SMB', 'HML', 'UMD']]
                X2 = tmpData[['Mkt-RF']]
                X3 = tmpData[['SP500-RF']]
                X1 = sm.add_constant(X1)
                X2 = sm.add_constant(X2)
                X3 = sm.add_constant(X3)
                model1 = sm.OLS(y, X1).fit()
                model2 = sm.OLS(y, X2).fit()
                model3 = sm.OLS(y, X3).fit()
                tmpDF1 = pd.DataFrame(model1.params).T
                tmpDF1.rename( columns = {'const' : 'alphaFF'}, inplace = True)
                tmpDF2 = pd.DataFrame(model2.params).T
                tmpDF2.rename( columns = {'const' : 'alphaCAPM', 'Mkt-RF' : 'Mkt-RF_only'}, inplace = True)
                tmpDF3 = pd.DataFrame(model3.params).T
                tmpDF3.rename( columns = {'const' : 'alphaSP500'}, inplace = True )
                tmpDF = pd.concat((tmpDF1, tmpDF2, tmpDF3), axis = 1)
                tmpDF['ticker'] = ticker
                Coef = Coef.append(tmpDF)
        Coef.to_csv(r'%s\Coef_%s_%s.csv' % (DataFolder, start_train.date(), end_train.date()), index = False)
        print 'Finished saving regression coefficients to: %s\Coef_%s_%s.csv' % (DataFolder, start_train.date(), end_train.date())
        return Coef

def regression(json_data, bandwidth):
    
    latency = []
    rtt_by_size = []

    # RTT object 
    # rtt = {[avg_rtt1]: [rtt1, rtt2, rtt3, ..., rtcx],
    #         [avg_rtt2]: [...]}

    for i in range(0, len(json_data)):
        latency.append(json_data[i]["latency"])
        rtt_by_size.append(json_data[i]["size"] * json_data[i]["rtt"])

    y = np.array(bandwidth).astype(np.float)
    z = np.array(latency).astype(np.float)
    r = np.array(rtt_by_size).astype(np.float)

    data = np.array([rtt_by_size, y])

    ones = np.ones(len(data[0]))
    X = sm.add_constant(np.column_stack((data[0], ones)))
    for ele in data[1:]:
        X = sm.add_constant(np.column_stack((ele, X)))
    results = sm.OLS(z, X).fit()
    print results.summary()

def reg_m(y, x):
    ones = np.ones(len(x[0]))
    X = sm.add_constant(np.column_stack((x[0], ones)))
    for ele in x[1:]:
        X = sm.add_constant(np.column_stack((ele, X)))
    results = sm.OLS(y, X).fit()
    return results

    def test_plot_influence(self, close_figures):
        infl = self.res.get_influence()
        fig = influence_plot(self.res)
        assert_equal(isinstance(fig, plt.Figure), True)
        # test that we have the correct criterion for sizes #3103
        try:
            sizes = fig.axes[0].get_children()[0]._sizes
            ex = sm.add_constant(infl.cooks_distance[0])
            ssr = sm.OLS(sizes, ex).fit().ssr
            assert_array_less(ssr, 1e-12)
        except AttributeError:
            import warnings
            warnings.warn('test not compatible with matplotlib version')

        fig = influence_plot(self.res, criterion='DFFITS')
        assert_equal(isinstance(fig, plt.Figure), True)
        try:
            sizes = fig.axes[0].get_children()[0]._sizes
            ex = sm.add_constant(np.abs(infl.dffits[0]))
            ssr = sm.OLS(sizes, ex).fit().ssr
            assert_array_less(ssr, 1e-12)
        except AttributeError:
            pass

        assert_raises(ValueError, influence_plot, self.res, criterion='unknown')

    def weak_instruments(self, n_sims=20):

        np.random.seed(1692)

        model = feedforward.FeedForwardModel(19, 1, dense_size=60, n_dense_layers=2)

        treatment_effects = []
        ols_betas, ols_ses = [], []
        old_corrs, new_corrs = [], []
        for _ in xrange(n_sims):
            df = self.treatment_gen.simulate_data(False)

            X = np.hstack((self.x, df['new_treat'].values[:, None]))
            Z = np.hstack((self.x, df['instrument'].values[:, None]))

            ols_beta, ols_se = self.fit_ols(df['treatment_effect'], X)
            ols_betas.append(ols_beta)
            ols_ses.append(ols_se)

            old_corr = df[['instrument', 'new_treat']].corr().values[0, 1]
            new_instrument, new_corr = model.fit_instruments(X, Z, df['treatment_effect'].values, batchsize=128)
            new_corrs.append(new_corr)
            old_corrs.append(old_corr)

            Z2 = Z.copy()
            Z2[:, -1] = new_instrument[:, 0]

            iv = IV2SLS(df['treatment_effect'].values.flatten(), add_constant(X), add_constant(Z2))

            model.reset_params()

        if new_corr:
            logger.info("Old corr: %.2f, New corr: %.2f", np.mean(old_corrs), np.mean(new_corrs))
        logger.info("Treatment effect (OLS): %.3f (%.4f)", np.mean(ols_betas), np.mean(ols_ses))
        logger.info("Treatment effect: %.3f (%.4f)", np.mean(treatment_effects), np.std(treatment_effects))

    def get_z_LinearRegression(self,xo,yo):
        print 'linear regression'
        dist_sigma = 1000
        xx= self.Knots[:,0]
        yy= self.Knots[:,1]
        dd = np.sqrt((xx-xo)**2 + (yy-yo)**2)
        print "dd",dd
        exponent = -(dd**2)/(2*(dist_sigma**2))
        print "exponent", exponent
        weights = np.exp(exponent)
        print "weights",weights

        X = self.Knots[:,0:2]
        X = sm.add_constant(X)
        y = self.Knots[:,2]

        mod_wls = sm.WLS(y, X, weights=weights)
        res_wls = mod_wls.fit()
        print(res_wls.summary())
        p = np.zeros((2,2),dtype=X.dtype)
        p[0,0] = xo
        p[0,1] = yo
        p[1,0] = xo
        p[1,0] = yo
        p = sm.add_constant(p)
        z = res_wls.predict(p)
        print "zshape",z

        return z[0]

def make_g_model(daily_results, daily_projections):
    daily_results_common = unify_dfs(daily_results)
    dfm = create_master(daily_results_common)
    dfm['NF'] = dfm['NF'].astype(float)
    dfm = eliminate_zeros(dfm)
    X = pd.get_dummies(dfm[['Salary', 'RG', 'NF', 'RW', 'POS', 'Depth']])
    X = pd.concat([X.drop('Depth',1), pd.get_dummies(X['Depth'])], 1)
    if 'POS_' in X.columns:
        X.drop('POS_', axis=1, inplace=True)
    #if 3 in X.columns:
        #X.drop(3, axis=1, inplace=True)
    print X.columns
    X = sm.add_constant(X)
    info = dfm[['Player', 'Date', 'Time']]
    y = dfm['FD']
    model=sm.OLS(y, X).fit()
    today = daily_projections[date_string]
    X = pd.get_dummies(today[['Salary', 'RG', 'NF', 'RW', 'POS', 'Depth']])
    X = pd.concat([X.drop('Depth',1), pd.get_dummies(X['Depth'])], 1)
    print X.columns
    if 'POS_' in X.columns:
        X.drop('POS_', axis=1, inplace=True)
    X = sm.add_constant(X)
    g_model = model.predict(X)
    return g_model  

def scatter(filename, x, y, line=True, xr=None, yr=None, x_title='', y_title='', title=None):
    if title is None:
        title = filename

    plt.figure(figsize=(24,18), dpi=600)
    plt.scatter(x, y)

    if xr is not None:
        plt.xlim(xr)
    if yr is not None:
        plt.ylim(yr)

    if line:
        est = sm.OLS(y, sm.add_constant(x)).fit()
        x_prime = np.linspace(min(x), max(x), 100)[:, np.newaxis]
        x_prime = sm.add_constant(x_prime)
        y_hat = est.predict(x_prime)
        line_plot1 = plt.plot(x_prime[:, 1], y_hat, 'r', alpha=0.9, label='r^2 = %s' % est.rsquared)
        #res = linregress(x,y)
        #line_plot2 = plt.plot([min(x), max(x)], [res[0]*min(x)+res[1], res[0]*max(x)+res[1]],
        #                      'g', alpha=0.9, label='r^2 = %s' % res[2])
        plt.legend(['r^2 = %s' % est.rsquared])

    plt.xlabel(x_title)
    plt.ylabel(y_title)
    plt.title(title)

    plt.savefig('%s.png' % filename, format='png')
    plt.savefig('%s.eps' % filename, format='eps')
    plt.close()

    def test_mvl_fuse_function(self):
        Y, D, P, T, G = generate_raw_samples()
        T = sm.add_constant(T, prepend=False)
        P = sm.add_constant(P, prepend=False)
        D = sm.add_constant(D, prepend=False)
        G = sm.add_constant(G, prepend=False)
        loo = LeaveOneOut(len(Y))
        er = []
        for train_idx, test_idx in loo:
            tm = taxi_view_model(train_idx, Y, T)
            pm = poi_view_model(train_idx, Y, P)
            gm = geo_view_model(train_idx, Y, G)
            dm = demo_view_model(train_idx, Y, D)
            models = [tm, pm, gm, dm]
            lm = mvl_fuse_function(models, train_idx, Y)
            
            
            tm_test = tm[0].predict(T[test_idx])
            pm_test = pm[0].predict(P[test_idx])
            gm_test = gm[0].predict(G[test_idx])
            dm_test = dm[0].predict(D[test_idx])
            
            newX_test = np.array([1, tm_test, pm_test, gm_test, dm_test])
            ybar = lm.predict(newX_test)
            y_error = ybar - Y[test_idx]
#            if np.abs(y_error / Y[test_idx]) > 0.8:
#                print test_idx, ybar, Y[test_idx], newX_test
            er.append(y_error)
        mre = np.mean(np.abs(er)) / np.mean(Y)
        print "MVL with linear fusion function MRE: {0}".format(mre)
        
        self.visualize_prediction_error(er, Y, "MVL linear combination")

def get_r_stat(sim_data):
    try:
        x = sm.add_constant(sim_data.rt_sampled.sort_values().reset_index().rt_sampled)
        y = sim_data.rt.sort_values().reset_index().rt
        fitted = sm.OLS(y,x).fit()
        log_x = sm.add_constant(sim_data.log_rt_sampled.sort_values().reset_index().log_rt_sampled)
        log_y = sim_data.log_rt.sort_values().reset_index().log_rt
        log_fitted = sm.OLS(log_y,log_x).fit()
        sub_out = pd.DataFrame([{'int_val': fitted.params.const,
                                'int_pval': fitted.pvalues.const,
                                'slope_val': fitted.params.rt_sampled,
                                'slope_pval':fitted.pvalues.rt_sampled,
                                'rsq': fitted.rsquared,
                                'rsq_adj': fitted.rsquared_adj,
                                'log_int_val': log_fitted.params.const,
                                'log_int_pval': log_fitted.pvalues.const,
                                'log_slope_val': log_fitted.params.log_rt_sampled,
                                'log_slope_pval': log_fitted.pvalues.log_rt_sampled,
                                'log_rsq': log_fitted.rsquared,
                                'log_rsq_adj': log_fitted.rsquared_adj}], index=[0])
    except:
        sub_out = pd.DataFrame([{'int_val': np.nan,
                                'int_pval': np.nan,
                                'slope_val': np.nan,
                                'slope_pval': np.nan,
                                'rsq': np.nan,
                                'rsq_adj': np.nan,
                                'log_int_val': np.nan,
                                'log_int_pval': np.nan,
                                'log_slope_val': np.nan,
                                'log_slope_pval': np.nan,
                                'log_rsq': np.nan,
                                'log_rsq_adj': np.nan}], index=[0])
    return sub_out

def models_pattern(data, matrix):
    y = np.array(data['survived'])
    ones = np.ones(len(matrix[0]))
    X = sm.add_constant(np.column_stack((matrix[0], ones)))
    for ele in matrix[1:]:
        X = sm.add_constant(np.column_stack((ele, X)))
    logit_model = sm.Logit(y, X)
    logit_res = logit_model.fit(maxiter=2000)
    print logit_res.summary()
    print logit_res.wald_test('1*x1 + 1*x2 + 1*x3')
    print
    probit_model = sm.Probit(y, X)
    probit_res = probit_model.fit(maxiter=2000)
    print probit_res.summary()
    print logit_res.wald_test('1*x1 + 1*x2 + 1*x3')
    print
    linear_model = sm.OLS(y, X)
    linear_res = linear_model.fit(maxiter=2000)
    result = 0.
    for array in X:
        for i, item in enumerate(array):
            result += linear_res.params[i] * item
    result /= (len(X))
    print 'Linear function value: {}'.format(result)
    print linear_res.summary()
    print linear_res.wald_test('1*x1 + 1*x2 + 1*x3')
    print

def detailedMultipleRegression(y, x):
    ones = np.ones(len(x[0]))
    X = sm.add_constant(np.column_stack((x[0], ones)))
    for ele in x[1:]:
        X = sm.add_constant(np.column_stack((ele, X)))
    results = sm.OLS(y, X).fit()
    return results # WARNING : coef not in right order !

def calcSScoreAgainstExisting(beta_mfac, beta_aux, df, n_aux):
  y = df.iloc[:, 0]
  X = df.iloc[:, 1:]
  resid_mfac = y - sm.add_constant(X).dot(beta_mfac)
  aux = calcAuxilaryArray(resid_mfac, n_aux)

  resid_aux = aux[1:] - sm.add_constant(aux[:-1]).dot(beta_aux)
  return calcSScore(beta_aux, resid_aux, aux[-1])

def get_recommend():

    # repeat process above to get data for recommendations
    data = flask.request.json
    headline = data["headline"]
    content = data["content"].encode('utf8')
    tags = data["tags"]
    day_published = data["day_pub"]
    channel = data["channel"]
    num_imgs = int(data["num_imgs"])
    index = [0]
    columns = [u'LDA_0_prob', u'LDA_1_prob', u'LDA_2_prob', u'LDA_3_prob',
       u'LDA_4_prob', u'LDA_5_prob', u'LDA_6_prob', u'LDA_7_prob',
       u'LDA_8_prob', u'LDA_9_prob', u'average_token_length_content',
       u'average_token_length_title', u'avg_negative_polarity',
       u'avg_positive_polarity', u'data_channel_is_bus',
       u'data_channel_is_entertainment', u'data_channel_is_lifestyle',
       u'data_channel_is_socmed', u'data_channel_is_tech',
       u'data_channel_is_world', u'global_grade_level',
       u'global_rate_negative_words', u'global_rate_positive_words',
       u'global_reading_ease', u'global_sentiment_abs_polarity',
       u'global_sentiment_polarity', u'global_subjectivity', u'is_weekend',
       u'max_abs_polarity', u'max_negative_polarity', u'max_positive_polarity',
       u'min_negative_polarity', u'min_positive_polarity', u'n_tokens_content',
       u'n_tokens_title', u'num_imgs', u'num_tags', u'num_videos',
       u'r_non_stop_unique_tokens', u'r_non_stop_words', u'r_unique_tokens',
       u'rate_negative_words', u'rate_positive_words',
       u'title_sentiment_abs_polarity', u'title_sentiment_polarity',
       u'title_subjectivity', u'weekday_is_friday',
       u'weekday_is_monday', u'weekday_is_saturday', u'weekday_is_sunday',
       u'weekday_is_thursday', u'weekday_is_tuesday', u'weekday_is_wednesday']
    data_df = pd.DataFrame(index=index, columns=columns)
    create_metadata_fields(data_df, num_imgs, tags, day_published, channel)
    create_NLP_features(data_df, headline, content)
    create_lda_features(data_df, content)
    results = {}
    create_metadata_fields(results, num_imgs, tags, day_published, channel)
    create_NLP_features(results, headline, content)
    create_lda_features(results, content)
    results['est_shares'] = round(pois_reg.predict(sm.add_constant(data_df))[0],-2)
    results['est_prob'] = round(RF_class.predict_proba(sm.add_constant(data_df))[0][1],2)

    # change day of week data to sunday
    data_df['weekday_is_monday'] = 0
    data_df['weekday_is_tuesday'] = 0
    data_df['weekday_is_wednesday'] = 0
    data_df['weekday_is_thursday'] = 0
    data_df['weekday_is_friday'] = 0
    data_df['weekday_is_saturday'] = 0
    data_df['weekday_is_sunday'] = 1
    data_df['is_weekend'] = 1

    # get results if change to sunday
    results['est_shares_sun'] = round(pois_reg.predict(sm.add_constant(data_df))[0],-2)
    results['est_prob_sun'] = round(RF_class.predict_proba(sm.add_constant(data_df))[0][1],2)

    # return results to javascript
    return flask.jsonify(results)