arma with exog

[joep]

is this a bug?

import numpy as np
import scikits.statsmodels.api as sm

exog = sm.add_constant(np.random.randn(100, 2), prepend=True)
endog = exog.sum(1) + 0.2 * np.random.randn(100)
modarma = sm.tsa.ARMA(endog, exog)
resarma = modarma.fit(order=(1,1))

Optimization terminated successfully.
         Current function value: 5.533772
         Iterations 12
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "sm_overview.py", line 36, in <module>
    resarma = modarma.fit(order=(1,1))
  File "E:\Josef\eclipsegworkspace\statsmodels-git\statsmodels-josef\scikits\statsmodels\tsa\arima_model.py", line 359, in fit
    start_params = self._fit_start_params((k_ar,k_ma,k))
  File "E:\Josef\eclipsegworkspace\statsmodels-git\statsmodels-josef\scikits\statsmodels\tsa\arima_model.py", line 76, in _fit_start_params
    start_params[:k] = ols_params
ValueError: shape mismatch: objects cannot be broadcast to a single shape

[Skipper Seabold (jsseabold)]

There is already a constant in exog, and then you try to add one with trend = 'c' in fit. The way, add_trend/add_constant are written, they will not add another one, but the model sets k_trend = 1. I guess the error message could be more informative. Maybe the helper functions should raise a warning or an error, "X already contains a constant"?

But I get a singular matrix error from the optimization if I try to fit it with trend='nc'. This could be the case when the model is badly mis-specified rather than a bug, but I'm not sure.

[Skipper Seabold (jsseabold)]

The docstring of add_trend says that it doesn't check for this, but it does check for an existing constant if trend=='c'. I guess this should be made more robust. A check for a trend would be something like

if np.any(np.all(np.diff(X,1,0)==1,0)), then has_trend.

[Skipper Seabold (jsseabold)]

Yeah, it just doesn't like the mis-specified model. Maybe we could catch these LinAlgErrors and warn that the model is probably mis-specified?

[~/]
[36]: resarma = modarma.fit(order=(1,0),trend='nc', disp=0)
RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N = 4 M = 12
 This problem is unconstrained.

           * * *

Tit = total number of iterations
Tnf = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip = number of BFGS updates skipped
Nact = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F = final function value

           * * *

   N Tit Tnf Tnint Skip Nact Projg F
    4 14 18 1 0 0 1.137D-05 -2.906D+01

CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH

 Total User time 1.500E-01 seconds.

[~/]
[37]: resarma.params
[37]: array([ 0.98864285, 1.00004556, 1.00576764, -0.0356923 ])

[joep (josef-pktd)]

I was just checking if I can run all models in the same or similar way

modols = sm.OLS(endog, exog)
modglsar = sm.GLSAR(endog, exog, rho=2)
modrlm = sm.RLM(endog, exog)
modlog = sm.Logit(endogbin, exog)
modglm = sm.GLM(endogbin, exog, family=sm.families.Binomial())
modar = sm.tsa.AR(endog)
modarma = sm.tsa.ARMA(endog)
modvar = sm.tsa.VAR(np.column_stack((endog, exog[:,1:])))

resols = modols.fit()
resglsar = modglsar.fit()
resrlm = modrlm.fit()
reslog = modlog.fit()
resglm = modglm.fit()
resar = modar.fit()
resarma = modarma.fit(order=(1,0))
resvar = modvar.fit()

ARMA is the only one without a default for everything in fit() but GenericLikelihood might also require something

I haven't looked at any details

instead of if np.any(np.all(np.diff(X,1,0)==1,0)), then has_trend.
maybe a check that the np.diff(X,1,0).var() < epsilon

for polynomial trends, I switched to rescaling trend to [-1, 1] or something like this, trend = np.linspace(-1,1,nobs)