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()