Estoy tratando de predecir una serie de tiempo en el paquete ARIMA de Python statsmodels con la inclusión de una variable exógena, pero no puedo encontrar la forma correcta de insertar la variable exógena en el paso de predicción. Veraquí para documentos

import numpy as np
from scipy import stats
import pandas as pd

import statsmodels.api as sm

vals = np.random.rand(13)
ts = pd.TimeSeries(vals)
df = pd.DataFrame(ts, columns=["test"])
df.index = pd.Index(pd.date_range("2011/01/01", periods = len(vals), freq = 'Q'))

fit1 = sm.tsa.ARIMA(df, (1,0,0)).fit()
#this works fine:
pred1 = fit1.predict(start=12, end = 16)
print(pred1)

Out[32]: 
2014-03-31    0.589121
2014-06-30    0.747575
2014-09-30    0.631322
2014-12-31    0.654858
2015-03-31    0.650093
Freq: Q-DEC, dtype: float64

ahora agregue una variable exógena de tendencia

exogx = np.array(range(1,14))
#to make this easy, let's look at the ols of the trend (arima(0,0,0))
fit2 = sm.tsa.ARIMA(df, (0,0,0),exog = exogx).fit()
print(fit2.params)

const    0.555226
x1       0.013132
dtype: float64

print(fit2.fittedvalues)

2011-03-31    0.568358
2011-06-30    0.581490
2011-09-30    0.594622
2011-12-31    0.607754
2012-03-31    0.620886
2012-06-30    0.634018
2012-09-30    0.647150
2012-12-31    0.660282
2013-03-31    0.673414
2013-06-30    0.686546
2013-09-30    0.699678
2013-12-31    0.712810
2014-03-31    0.725942
Freq: Q-DEC, dtype: float64

Tenga en cuenta que, como es de esperar, esta es una línea de tendencia, que aumenta 0.013132 con cada incremento en el tiempo (por supuesto, estos son datos aleatorios, por lo que si lo ejecuta, los valores serán diferentes, pero la historia de tendencia positiva o negativa será la mismo). Entonces, el siguiente valor (para tiempo = 14) debe ser 0.555226 + 0.013132 * 14 = 0.739074.

#out of sample exog should be (14,15,16)
pred2 = fit2.predict(start = 12, end = 16, exog = np.array(range(13,17)))
print(pred2)
2014-03-31    0.725942
2014-06-30    0.568358
2014-09-30    0.581490
2014-12-31    0.594622
2015-03-31    0.765338
Freq: Q-DEC, dtype: float64

Entonces, el 31-03-2014 predice (la última muestra) correctamente, pero el 30/06/2014 comienza de nuevo desde el principio (t = 1), pero observe el 31-03-2015 (en realidad, siempre la última observación del pronóstico, independientemente del horizonte) recoge t = 16 (es decir, (valor - intercepción) / beta = (0.765338 - 0.555226) /0.013132).

Para aclarar esto, observe lo que sucede cuando inflo los valores de la estera x

fit2.predict(start = 12, end = 16, exog = np.array(range(13,17))*10000)
Out[41]: 
2014-03-31       0.725942
2014-06-30       0.568358
2014-09-30       0.581490
2014-12-31       0.594622
2015-03-31    2101.680532
Freq: Q-DEC, dtype: float64

¿Ves que el 31-03-2015 explotó, pero ninguno de los otros valores xmat fueron considerados? ¿¿¿Qué estoy haciendo mal aquí???

He intentado jugar con todas las formas en que sé cómo pasar la variable exog (cambiar la dimensión, hacer que el exog sea una matriz, hacer que el exog sea tan largo como la entrada más el horizonte, etc., etc.). Cualquier sugerencia sería muy apreciada.

Estoy usando 2.7 de Anaconda2.1 numpy 1.8.1 scipy 0.14.0 pandas 0.14.0 statsmodels 0.5.0

y he verificado el problema en Windows 7 de 64 bits y centos de 64 bits.

Además, algunas cosas. Estoy usando ARIMA para la funcionalidad de ARIMA y lo anterior es solo para ilustración (es decir, no puedo "simplemente usar OLS ...", como imagino que se sugerirá). Tampoco puedo "simplemente usar R" debido a las restricciones del proyecto (y, en general, la falta de soporte de R en la base de Spark).

Aquí están todas las partes interesantes del código en caso de que quiera probarlo usted mismo.

import numpy as np
from scipy import stats
import pandas as pd
import statsmodels.api as sm

vals = np.random.rand(13)
ts = pd.TimeSeries(vals)
df = pd.DataFrame(ts, columns=["test"])
df.index = pd.Index(pd.date_range("2011/01/01", periods = len(vals), freq = 'Q'))

exogx = np.array(range(1,14))
fit2 = sm.tsa.ARIMA(df, (0,0,0),exog = exogx).fit()
print(fit2.fittedvalues)
pred2 = fit2.predict(start = 12, end = 16, exog = np.array(range(13,17))*10000)
print(pred2)