Variável exógena Python ARIMA fora da amostra

Estou tentando prever uma série temporal no pacote ARIMA do python statsmodels com a inclusão de uma variável exógena, mas não consigo descobrir a maneira correta de inserir a variável exógena na etapa de previsão. Vejoaqui 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

agora adicione uma variável exógena de tendência

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

Observe que, como seria de esperar, esta é uma linha de tendência, aumentando 0,013132 a cada aumento de escala no tempo (é claro que são dados aleatórios, portanto, se você executá-lo, os valores serão diferentes, mas a história de tendência positiva ou negativa será a mesmo). Portanto, o próximo valor (para tempo = 14) deve 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

Portanto, 31/03/2014 prevê (a última amostra) corretamente, mas 30/06/2014 começa no início (t = 1), mas observe 31/03/2015 (na verdade, sempre a última observação da previsão, independentemente do horizonte) pega t = 16 (ou seja, (valor - interceptação) / beta = (0,765338 - 0,555226) / 0,013132).

Para deixar isso mais claro, observe o que acontece quando eu inflo os valores de x mat

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

Veja que 31/03/2015 explodiu, mas nenhum dos outros valores de xmat foram considerados? O que eu estou fazendo errado aqui???

Eu tentei brincar de todas as formas que sei para passar na variável exog (mudar de dimensão, transformar o exog em uma matriz, fazer o exog contanto que a entrada seja mais o horizonte, etc, etc, etc). Qualquer sugestão seria muito apreciada.

Estou usando o 2.7 do Anaconda2.1 numpy 1.8.1 scipy 0.14.0 pandas 0.14.0 statsmodels 0.5.0

e verificamos o problema no windows 7 de 64 bits e no centos de 64 bits.

Além disso, algumas coisas. Estou usando o ARIMA para a funcionalidade ARIMA e o exemplo acima é apenas para ilustração (ou seja, não posso "apenas usar OLS ...", como imagino que será sugerido). Também não posso "apenas usar R" devido às restrições do projeto (e mais geralmente, à falta de suporte do R no Spark base).

Aqui estão todas as partes interessantes do código, caso você queira tentar por conta própria

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)

questionAnswers(3)

yourAnswerToTheQuestion