def getSeries(sData):
    baseUrl = "http://analisi.ad.mediamond.it/jsonp.php?"
    resq = requests.get(baseUrl+urllib.urlencode(sData))##,'callback':'?'}))
    ser = np.array(resq.json()['data']).astype(np.float)
    timeL = pd.DataFrame({'t':(np.array(resq.json()['data'])[:,0])/1000-3600})
    timeL = timeL.t.apply(lambda x: datetime.datetime.fromtimestamp(x))
    trainD = pd.DataFrame({'y':ser[:,1]/1000000},index=[timeL[0]+datetime.timedelta(x) for x in range(timeL.size)])
    trainD['t'] = ser[:,0]/1000000000
    trainD.dropna(inplace=True)
    tempD = trainD.copy(deep=True)
    tempD['week'] = [x.isocalendar()[1] for x in tempD.index]
    sWeek = tempD.groupby(["week"]).sum()
    sWeek['count'] = tempD.groupby(["week"]).size()
    sWeek['y'] = sWeek['y']/sWeek['count']
    normY = sWeek['y'].max()
    sWeek['y'] = sWeek['y']/normY
    base = datetime.datetime.strptime("2017-01-01","%Y-%m-%d")
    nWeek = sWeek.shape[0]*7
    sWeek.index = [base + datetime.timedelta(days=x) for x in range(0,nWeek,7)]
    sWeek.t = [float(calendar.timegm(x.utctimetuple()))/1000000 for x in sWeek.index]
    sWeek['month'] = sWeek.index.month
    sMonth = sWeek.groupby(['month']).sum()
    sMonth['count'] = sWeek.groupby(['month']).size()
    sMonth['y'] = sMonth['y']/sMonth['count']
    normY = sMonth['y'].max()
    sMonth['y'] = sMonth['y']/normY
    sMonth.index = ['2017-%02d-01' % x for x in range(1,len(sMonth)+1)]
    sMonth.index = [datetime.datetime.strptime(x,"%Y-%m-%d") for x in sMonth.index]
    sMonth.t = [float(calendar.timegm(x.utctimetuple()))/1000000 for x in sMonth.index]
    trainD['roll'] = smooth(trainD['y'],4,11)
    trainD['e_av'] = pd.ewma(trainD['y'],halflife=14)
    #trainD['e_av'] = pd.Series.ewm(trainD['y'],halflife=14)
    trainD['dev'] = trainD['y'] - trainD['y'].shift()
    trainD['diff'] = trainD['y'] - trainD['e_av']
    trainD['stat'] = trainD['dev']
    return trainD, sWeek, sMonth
			

Fetch series


def extSeries(sDay,nAhead,x0,hWeek):
    nLin = sDay.shape[0] + nAhead
    nFit = sDay.shape[0] if int(x0[6]) <= 14 else int(x0[6])
    sDay['hist'] = sp.interpolate.interp1d(hWeek.t,hWeek.y,kind="cubic")(sDay['t'])
    sDay['hist'] = sDay['hist']/sDay['hist'].mean()
    lmFor = 'e_av ~ 1 + t + I(t**2) + I(t**3) + I(t**4) + I(t**5)'
    lm = smf.ols(formula=lmFor,data=sDay.tail(nFit)).fit()
    sDay['stat'] = (sDay['y']*sDay['hist']*x0[10]-lm.predict(sDay))
    #print 2.*np.pi/(sDay.t[1]-sDay.t[0])
    t_test = np.linspace(sDay['t'][0],sDay['t'][sDay.shape[0]-1]+sDay.t[nAhead]-sDay.t[0],nLin)
    testD = pd.DataFrame({'t':t_test},index=[sDay.index[0]+datetime.timedelta(days=x) for x in range(nLin)])
    testD['y'] = sDay.y
    testD['hist'] = sp.interpolate.interp1d(hWeek.t,hWeek.y,kind="cubic")(testD['t'])
    testD['hist'] = testD['hist']/sDay['hist'].mean()
    testD['trend'] = lm.predict(testD)
    testD = testD.tail(nFit+nAhead)
    freqP = [x0[7],x0[8]]
    #print 2.*np.pi/(sDay.t[7]-sDay.t[0])
    def fun(x,t):
        return x[0] + x[1] * np.sin(freqP[0]*t + x[2])*(1 + x[3]*np.sin(freqP[1]*t + x[4]))    ##confInt = stats.t.interval(0.95,len(y)-1,loc=np.mean(y),scale=stats.sem(y))
    def fun_min(x,t,y):
        return fun(x,t) - y
    res_lsq = least_squares(fun_min,x0,args=(sDay.t,sDay.stat))
    ##res_robust = least_squares(fun_min,x0,loss='soft_l1',f_scale=0.1,args=(sDay.t,sDay.stat))
    testD['lsq'] = fun(res_lsq[0],testD.t) # fun(res_robust.x,t_test)
    sDay['resid'] = sDay['y'] - fun(res_lsq[0],sDay.t)/x0[10] - lm.predict(sDay)
    rSquare = (sDay['resid'] - sDay['resid'].mean())**2
    testD['pred'] = (testD['lsq']/x0[10] + lm.predict(testD))
    testD = testD.drop(testD.index[0])
    return testD, res_lsq[0], rSquare.sum()
			

Least squares


def serArma(sDay,nAhead,x0,hWeek):
    x1 = [int(x) for x in x0]
    sDay['hist'] = sp.interpolate.interp1d(hWeek.t,hWeek.y,kind="cubic")(sDay['t'])
    lmFor = 'e_av ~ 1 + t + I(t**2) + I(t**3) + I(t**4) + I(t**5)'
    lm = smf.ols(formula=lmFor,data=sDay).fit()
    sDay['hist'] = sDay['hist']/sDay['hist'].mean()
    sDay['stat'] = (sDay['y']*sDay['hist']*x0[10]-lm.predict(sDay))
    dta = sDay['y']
    dta.index = [pd.datetime.strptime(str(x)[0:10],'%Y-%m-%d') for x in dta.index]
    t_line = [float(calendar.timegm(x.utctimetuple()))/1000000 for x in dta.index]
    t_line1 = [float(calendar.timegm(x.utctimetuple()))/1000000 for x in hWeek.index]
    sExog = pd.DataFrame({'y':sp.interpolate.interp1d(t_line1,hWeek.y,kind="cubic")(t_line)})
    #par = bestArima(dta,sExog)
    sExog.index = dta.index
    result = sm.tsa.ARIMA(dta,(x1[0],x1[1],x1[2])).fit(trend="c",method='css-mle',exog=sExog)
    predT = [str(dta.index[0])[0:10],str(dta.index[len(dta)-1])[0:10]]
    histS = pd.DataFrame({'pred':result.predict(start=predT[0],end=predT[1])})
    sDay['resid'] = sDay['y'] - histS['pred']
    rSquare = (sDay['resid'] - sDay['resid'].mean())**2
    predT = [str(dta.index[0])[0:10],str(dta.index[len(dta)-1]+datetime.timedelta(days=nAhead))[0:10]]    
    predS = pd.DataFrame({'pred':result.predict(start=predT[0],end=predT[1])})
    predS.index = [dta.index[0]+datetime.timedelta(days=x) for x in range(0,len(dta)+nAhead)]
    predS['t'] = [float(calendar.timegm(x.utctimetuple()))/1000000 for x in predS.index]
    predS['hist'] = sp.interpolate.interp1d(t_line1,hWeek.y,kind="cubic")(predS['t'])
    predS['hist'] = predS['hist']/predS['hist'].mean()
    predS['pred'] = predS['pred']*predS['hist']*x0[10]
    predS['trend'] = lm.predict(predS)
    predS['lsq'] = 0
    ##predS['pred'] = (predS['pred']*predS['hist']+predS['trend'])
    x2 = x0#.values
    x2[3:(3+len(result.params))] = result.params
    x2[9] = 1
    predS = predS.drop(predS.index[0])
    return predS, pd.DataFrame(x2), rSquare.sum()
			

Arma forecast


def serHolt(sDay,nAhead,x0,hWeek):
    Y = [x for x in sDay.y]
    ##Yht, alpha, beta, gamma, rmse = ht.additive([x for x in sDay.y],int(x0[0]),nAhead,x0[1],x0[2],x0[3])
    nAv = int(x0[0]) if int(x0[0]) > 1 else 5
    Yht, alpha, beta, gamma, rmse = ht.additive([x for x in sDay.y],nAv,nAhead,x0[1],x0[2],x0[3])
    sDay['resid'] = sDay['y'] - Yht[0:sDay.shape[0]]
    x0[1] = alpha
    x0[2] = beta
    x0[3] = gamma
    x0[4] = rmse
    nLin = sDay.shape[0] + nAhead
    t_test = np.linspace(sDay['t'][0],sDay['t'][sDay.shape[0]-1]+sDay.t[nAhead]-sDay.t[0],nLin)
    predS = pd.DataFrame({'t':t_test},index=[sDay.index[0]+datetime.timedelta(days=x) for x in range(nLin)])
    predS['pred'] = Yht
    predS['hist'] = sp.interpolate.interp1d(hWeek.t,hWeek.y,kind="cubic")(predS['t'])
    predS['hist'] = predS['hist']/predS['hist'].mean()
    predS['pred'] = predS['pred']*predS['hist']*x0[10]
    lmFor = 'e_av ~ 1 + t + I(t**2) + I(t**3) + I(t**4) + I(t**5)'
    lm = smf.ols(formula=lmFor,data=sDay).fit()
    predS['trend'] = lm.predict(predS)*x0[4]
    predS['lsq'] = 0
    sDay['resid'] = sDay['y'] - predS['pred'][0:sDay.shape[0]]
    rSquare = (sDay['resid'] - sDay['resid'].mean())**2
    return predS, x0, rSquare.sum()
			

Holt winters


def serAuto(sDay,nAhead,x0,hWeek):
    sDay['hist'] = sp.interpolate.interp1d(hWeek.t,hWeek.y,kind="cubic")(sDay['t'])
    lmFor = 'e_av ~ 1 + t + I(t**2) + I(t**3) + I(t**4) + I(t**5)'
    lm = smf.ols(formula=lmFor,data=sDay).fit()
    sDay['stat'] = (sDay['y']*sDay['hist']*x0[10]-lm.predict(sDay))
    sDay['hist'] = sDay['hist']/sDay['hist'].mean()
    todayD = datetime.datetime.today()
    todayD = todayD.replace(hour=0,minute=0,second=0,microsecond=0)
    dta = pd.DataFrame({'y':sDay.y})
    dta['day'] = sDay.index.weekday
    phase = dta.head(int(x0[6])).groupby(['day']).mean()
    phase['std'] = dta.groupby(['day']).std()['y']
    phase = phase.sort_values(['y'],ascending=False)
    phase['csum'] = phase['y'].cumsum()/phase['y'].sum()
    phaseN = phase.index[0] - todayD.weekday()
    r,q,p = sm.tsa.acf(sDay['y'].tail(phaseN+int(x0[6])).squeeze(),qstat=True)
    def fit_fun(x,decay):
        return np.exp(-decay*x)
    popt, pcov = curve_fit(fit_fun,np.array(range(0,6)),r[0:6]-min(r),p0=(x0[5]))
    X = np.array(range(0,r.size,7))
    popt1, pcov1 = curve_fit(fit_fun,X,r[X],p0=(x0[5]))
    autD = pd.DataFrame({'r':r,'exp':fit_fun(range(0,r.size),popt),'exp1':fit_fun(range(0,r.size),popt1)})    
    x0[5] = popt
    testD = pd.DataFrame(index=[todayD + datetime.timedelta(days=x) for x in range(-sDay.shape[0],nAhead)])
    testD['t'] = [float(calendar.timegm(x.utctimetuple()))/1000000. for x in testD.index]
    wN = 0
    sY = np.random.normal(phase['y'].head(1),dta.y.std())
    testD['hist'] = sp.interpolate.interp1d(hWeek.t,hWeek.y,kind="cubic")(testD['t'])
    testD['pred'] = 0
    for i in testD.index:
        wN = 6 - np.abs(phase.index[0] - i.weekday())
        wN = wN + 1 if wN < 6 else 0
        if(wN == 0):
            sY = np.random.normal(phase['y'].head(1),dta.y.std()/2)
        sY = sY*(1+testD['hist'][i]*x0[10])
        testD.loc[i,'pred'] = sY*fit_fun(float(wN),popt)

    testD['pred1'] = testD['pred']
    testD['pred'] = smooth(testD['pred'],16,5)
    sDay['resid'] = sDay['y'] - testD['pred'][0:sDay.shape[0]]
    sDay['resid1'] = sDay['resid']
    freqP = [x0[7],x0[8]]
    def fun(x,t):
        return x[0] + x[1] * np.sin(freqP[0]*t + x[2])*(1 + x[3]*np.sin(freqP[1]*t + x[4]))    
    def fun_min(x,t,y):
        return fun(x,t) - y
    res_lsq = least_squares(fun_min,x0,args=(sDay['t'],sDay['resid']))
    testD['lsq'] = fun(res_lsq[0],testD['t']) # fun(res_robust.x,t_test)
    x0[0:res_lsq[0].size] = res_lsq[0]
    testD['pred2'] = testD['pred']
    testD['pred'] = testD['pred'] + testD['lsq']
    sDay['resid'] = sDay['y'] - testD['pred'][0:sDay.shape[0]]    
    rSquare = (sDay['resid'] - sDay['resid'].mean())**2
    testD['trend'] = lm.predict(testD)
    # sDay.to_csv('tmpAuto1.csv')
    # testD.to_csv('tmpAuto2.csv')
    # autD.to_csv('tmpAuto3.csv')
    return testD, x0, rSquare.sum()

			

Auto correlated simulations


def SerBayes(sDay,nAhead,x0,hWeek):
    x1 = [int(x) for x in x0]
    dta = sDay['y']
    dta.index = [pd.datetime.strptime(str(x)[0:10],'%Y-%m-%d') for x in dta.index]
    t_line = [float(calendar.timegm(x.utctimetuple()))/1000000 for x in dta.index]
    dta.index = t_line
    model = pydlm.dlm(dta)
    model = model + pydlm.trend(degree=1,discount=0.98,name='a',w=10.0)
    model = model + pydlm.dynamic(features=[[v] for v in t_line],discount=1,name='b',w=10.0)
    model = model + pydlm.autoReg(degree=3,data=dta.values,name='ar3',w=1.0)
    allStates = model.getLatentState(filterType='forwardFilter')
    model.evolveMode('independent')
    model.noisePrior(2.0)
    model.fit()
    model.plot()
    model.turnOff('predict')
    model.plotCoef(name='a')
    model.plotCoef(name='b')
    model.plotCoef(name='ar3')
			

Bayesian