PageRenderTime 60ms CodeModel.GetById 18ms RepoModel.GetById 0ms app.codeStats 1ms

/statsmodels/tsa/tests/test_arima.py

http://github.com/statsmodels/statsmodels
Python | 2746 lines | 2653 code | 79 blank | 14 comment | 4 complexity | 4e2322760d7505b95873cdee2fe1766c MD5 | raw file
Possible License(s): BSD-3-Clause
  1. from statsmodels.compat.platform import (PLATFORM_OSX, PLATFORM_WIN,
  2. PLATFORM_WIN32)
  3. from statsmodels.compat.python import lrange
  4. import os
  5. import pickle
  6. import warnings
  7. from io import BytesIO
  8. import numpy as np
  9. import pandas as pd
  10. import pytest
  11. from numpy.testing import assert_almost_equal, assert_allclose, assert_raises
  12. from pandas import DatetimeIndex, date_range, period_range
  13. import statsmodels.sandbox.tsa.fftarma as fa
  14. from statsmodels.datasets.macrodata import load_pandas as load_macrodata_pandas
  15. from statsmodels.regression.linear_model import OLS
  16. from statsmodels.tools.sm_exceptions import (
  17. ValueWarning, HessianInversionWarning, SpecificationWarning,
  18. MissingDataError)
  19. from statsmodels.tools.testing import assert_equal
  20. from statsmodels.tsa.ar_model import AutoReg
  21. from statsmodels.tsa.arima_model import ARMA, ARIMA
  22. from statsmodels.tsa.arima_process import arma_generate_sample
  23. from statsmodels.tsa.arma_mle import Arma
  24. from statsmodels.tsa.tests.results import results_arma, results_arima
  25. DECIMAL_4 = 4
  26. DECIMAL_3 = 3
  27. DECIMAL_2 = 2
  28. DECIMAL_1 = 1
  29. current_path = os.path.dirname(os.path.abspath(__file__))
  30. ydata_path = os.path.join(current_path, 'results', 'y_arma_data.csv')
  31. with open(ydata_path, "rb") as fd:
  32. y_arma = np.genfromtxt(fd, delimiter=",", skip_header=1, dtype=float)
  33. cpi_dates = period_range(start='1959q1', end='2009q3', freq='Q')
  34. sun_dates = period_range(start='1700', end='2008', freq='A')
  35. cpi_predict_dates = period_range(start='2009q3', end='2015q4', freq='Q')
  36. sun_predict_dates = period_range(start='2008', end='2033', freq='A')
  37. def test_compare_arma():
  38. # this is a preliminary test to compare arma_kf, arma_cond_ls
  39. # and arma_cond_mle
  40. # the results returned by the fit methods are incomplete
  41. # for now without random.seed
  42. np.random.seed(9876565)
  43. x = fa.ArmaFft([1, -0.5], [1., 0.4], 40).generate_sample(nsample=200,
  44. burnin=1000)
  45. modkf = ARMA(x, (1, 1))
  46. reskf = modkf.fit(trend='nc', disp=-1)
  47. dres = reskf
  48. modc = Arma(x)
  49. resls = modc.fit(order=(1, 1))
  50. rescm = modc.fit_mle(order=(1, 1), start_params=[0.4, 0.4, 1.], disp=0)
  51. # decimal 1 corresponds to threshold of 5% difference
  52. # still different sign corrcted
  53. assert_almost_equal(resls[0] / dres.params, np.ones(dres.params.shape),
  54. decimal=1)
  55. # rescm also contains variance estimate as last element of params
  56. assert_almost_equal(rescm.params[:-1] / dres.params,
  57. np.ones(dres.params.shape), decimal=1)
  58. class CheckArmaResultsMixin(object):
  59. """
  60. res2 are the results from gretl. They are in results/results_arma.
  61. res1 are from statsmodels
  62. """
  63. decimal_params = DECIMAL_4
  64. def test_params(self):
  65. assert_almost_equal(self.res1.params, self.res2.params,
  66. self.decimal_params)
  67. decimal_aic = DECIMAL_4
  68. def test_aic(self):
  69. assert_almost_equal(self.res1.aic, self.res2.aic, self.decimal_aic)
  70. decimal_bic = DECIMAL_4
  71. def test_bic(self):
  72. assert_almost_equal(self.res1.bic, self.res2.bic, self.decimal_bic)
  73. decimal_arroots = DECIMAL_4
  74. def test_arroots(self):
  75. assert_almost_equal(self.res1.arroots, self.res2.arroots,
  76. self.decimal_arroots)
  77. decimal_maroots = DECIMAL_4
  78. def test_maroots(self):
  79. assert_almost_equal(self.res1.maroots, self.res2.maroots,
  80. self.decimal_maroots)
  81. decimal_bse = DECIMAL_2
  82. def test_bse(self):
  83. assert_almost_equal(self.res1.bse, self.res2.bse, self.decimal_bse)
  84. decimal_cov_params = DECIMAL_4
  85. def test_covparams(self):
  86. assert_almost_equal(self.res1.cov_params(), self.res2.cov_params,
  87. self.decimal_cov_params)
  88. decimal_hqic = DECIMAL_4
  89. def test_hqic(self):
  90. assert_almost_equal(self.res1.hqic, self.res2.hqic, self.decimal_hqic)
  91. decimal_llf = DECIMAL_4
  92. def test_llf(self):
  93. assert_almost_equal(self.res1.llf, self.res2.llf, self.decimal_llf)
  94. decimal_resid = DECIMAL_4
  95. def test_resid(self):
  96. assert_almost_equal(self.res1.resid, self.res2.resid,
  97. self.decimal_resid)
  98. decimal_fittedvalues = DECIMAL_4
  99. def test_fittedvalues(self):
  100. assert_almost_equal(self.res1.fittedvalues, self.res2.fittedvalues,
  101. self.decimal_fittedvalues)
  102. decimal_pvalues = DECIMAL_2
  103. def test_pvalues(self):
  104. assert_almost_equal(self.res1.pvalues, self.res2.pvalues,
  105. self.decimal_pvalues)
  106. decimal_t = DECIMAL_2 # only 2 decimal places in gretl output
  107. def test_tvalues(self):
  108. assert_almost_equal(self.res1.tvalues, self.res2.tvalues,
  109. self.decimal_t)
  110. decimal_sigma2 = DECIMAL_4
  111. def test_sigma2(self):
  112. assert_almost_equal(self.res1.sigma2, self.res2.sigma2,
  113. self.decimal_sigma2)
  114. @pytest.mark.smoke
  115. def test_summary(self):
  116. self.res1.summary()
  117. @pytest.mark.smoke
  118. def test_summary2(self):
  119. self.res1.summary2()
  120. class CheckForecastMixin(object):
  121. decimal_forecast = DECIMAL_4
  122. def test_forecast(self):
  123. assert_almost_equal(self.res1.forecast_res, self.res2.forecast,
  124. self.decimal_forecast)
  125. decimal_forecasterr = DECIMAL_4
  126. def test_forecasterr(self):
  127. assert_almost_equal(self.res1.forecast_err, self.res2.forecasterr,
  128. self.decimal_forecasterr)
  129. class CheckDynamicForecastMixin(object):
  130. decimal_forecast_dyn = 4
  131. def test_dynamic_forecast(self):
  132. assert_almost_equal(self.res1.forecast_res_dyn, self.res2.forecast_dyn,
  133. self.decimal_forecast_dyn)
  134. def test_forecasterr(self):
  135. assert_almost_equal(self.res1.forecast_err_dyn,
  136. self.res2.forecasterr_dyn,
  137. DECIMAL_4)
  138. class CheckArimaResultsMixin(CheckArmaResultsMixin):
  139. def test_order(self):
  140. assert self.res1.k_diff == self.res2.k_diff
  141. assert self.res1.k_ar == self.res2.k_ar
  142. assert self.res1.k_ma == self.res2.k_ma
  143. decimal_predict_levels = DECIMAL_4
  144. def test_predict_levels(self):
  145. assert_almost_equal(self.res1.predict(typ='levels'), self.res2.linear,
  146. self.decimal_predict_levels)
  147. class Test_Y_ARMA11_NoConst(CheckArmaResultsMixin, CheckForecastMixin):
  148. @classmethod
  149. def setup_class(cls):
  150. endog = y_arma[:, 0]
  151. cls.res1 = ARMA(endog, order=(1, 1)).fit(trend='nc', disp=-1)
  152. (cls.res1.forecast_res, cls.res1.forecast_err,
  153. confint) = cls.res1.forecast(10)
  154. cls.res2 = results_arma.Y_arma11()
  155. def test_pickle(self):
  156. fh = BytesIO()
  157. # test wrapped results load save pickle
  158. self.res1.save(fh)
  159. fh.seek(0, 0)
  160. res_unpickled = self.res1.__class__.load(fh)
  161. assert type(res_unpickled) is type(self.res1) # noqa: E721
  162. class Test_Y_ARMA14_NoConst(CheckArmaResultsMixin):
  163. @classmethod
  164. def setup_class(cls):
  165. endog = y_arma[:, 1]
  166. cls.res1 = ARMA(endog, order=(1, 4)).fit(trend='nc', disp=-1)
  167. cls.res2 = results_arma.Y_arma14()
  168. @pytest.mark.slow
  169. class Test_Y_ARMA41_NoConst(CheckArmaResultsMixin, CheckForecastMixin):
  170. @classmethod
  171. def setup_class(cls):
  172. endog = y_arma[:, 2]
  173. cls.res1 = ARMA(endog, order=(4, 1)).fit(trend='nc', disp=-1)
  174. (cls.res1.forecast_res, cls.res1.forecast_err,
  175. confint) = cls.res1.forecast(10)
  176. cls.res2 = results_arma.Y_arma41()
  177. cls.decimal_maroots = DECIMAL_3
  178. class Test_Y_ARMA22_NoConst(CheckArmaResultsMixin):
  179. @classmethod
  180. def setup_class(cls):
  181. endog = y_arma[:, 3]
  182. cls.res1 = ARMA(endog, order=(2, 2)).fit(trend='nc', disp=-1)
  183. cls.res2 = results_arma.Y_arma22()
  184. class Test_Y_ARMA50_NoConst(CheckArmaResultsMixin, CheckForecastMixin):
  185. @classmethod
  186. def setup_class(cls):
  187. endog = y_arma[:, 4]
  188. cls.res1 = ARMA(endog, order=(5, 0)).fit(trend='nc', disp=-1)
  189. (cls.res1.forecast_res, cls.res1.forecast_err,
  190. confint) = cls.res1.forecast(10)
  191. cls.res2 = results_arma.Y_arma50()
  192. class Test_Y_ARMA02_NoConst(CheckArmaResultsMixin):
  193. @classmethod
  194. def setup_class(cls):
  195. endog = y_arma[:, 5]
  196. cls.res1 = ARMA(endog, order=(0, 2)).fit(trend='nc', disp=-1)
  197. cls.res2 = results_arma.Y_arma02()
  198. class Test_Y_ARMA11_Const(CheckArmaResultsMixin, CheckForecastMixin):
  199. @classmethod
  200. def setup_class(cls):
  201. endog = y_arma[:, 6]
  202. cls.res1 = ARMA(endog, order=(1, 1)).fit(trend="c", disp=-1)
  203. (cls.res1.forecast_res, cls.res1.forecast_err,
  204. confint) = cls.res1.forecast(10)
  205. cls.res2 = results_arma.Y_arma11c()
  206. class Test_Y_ARMA14_Const(CheckArmaResultsMixin):
  207. @classmethod
  208. def setup_class(cls):
  209. endog = y_arma[:, 7]
  210. cls.res1 = ARMA(endog, order=(1, 4)).fit(trend="c", disp=-1)
  211. cls.res2 = results_arma.Y_arma14c()
  212. class Test_Y_ARMA41_Const(CheckArmaResultsMixin, CheckForecastMixin):
  213. @classmethod
  214. def setup_class(cls):
  215. endog = y_arma[:, 8]
  216. cls.res2 = results_arma.Y_arma41c()
  217. cls.res1 = ARMA(endog, order=(4, 1)).fit(trend="c", disp=-1,
  218. start_params=cls.res2.params)
  219. (cls.res1.forecast_res, cls.res1.forecast_err,
  220. confint) = cls.res1.forecast(10)
  221. cls.decimal_cov_params = DECIMAL_3
  222. cls.decimal_fittedvalues = DECIMAL_3
  223. cls.decimal_resid = DECIMAL_3
  224. cls.decimal_params = DECIMAL_3
  225. class Test_Y_ARMA22_Const(CheckArmaResultsMixin):
  226. @classmethod
  227. def setup_class(cls):
  228. endog = y_arma[:, 9]
  229. cls.res1 = ARMA(endog, order=(2, 2)).fit(trend="c", disp=-1)
  230. cls.res2 = results_arma.Y_arma22c()
  231. def test_summary(self):
  232. # regression test for html of roots table #4434
  233. # we ignore whitespace in the assert
  234. summ = self.res1.summary()
  235. summ_roots = """\
  236. <tableclass="simpletable">
  237. <caption>Roots</caption>
  238. <tr>
  239. <td></td><th>Real</th><th>Imaginary</th><th>Modulus</th><th>Frequency</th>
  240. </tr>
  241. <tr>
  242. <th>AR.1</th><td>1.0991</td><td>-1.2571j</td><td>1.6698</td><td>-0.1357</td>
  243. </tr>
  244. <tr>
  245. <th>AR.2</th><td>1.0991</td><td>+1.2571j</td><td>1.6698</td><td>0.1357</td>
  246. </tr>
  247. <tr>
  248. <th>MA.1</th><td>-1.1702</td><td>+0.0000j</td><td>1.1702</td><td>0.5000</td>
  249. </tr>
  250. <tr>
  251. <th>MA.2</th><td>1.2215</td><td>+0.0000j</td><td>1.2215</td><td>0.0000</td>
  252. </tr>
  253. </table>"""
  254. assert_equal(summ.tables[2]._repr_html_().replace(' ', ''),
  255. summ_roots.replace(' ', ''))
  256. class Test_Y_ARMA50_Const(CheckArmaResultsMixin, CheckForecastMixin):
  257. @classmethod
  258. def setup_class(cls):
  259. endog = y_arma[:, 10]
  260. cls.res1 = ARMA(endog, order=(5, 0)).fit(trend="c", disp=-1)
  261. (cls.res1.forecast_res, cls.res1.forecast_err,
  262. confint) = cls.res1.forecast(10)
  263. cls.res2 = results_arma.Y_arma50c()
  264. class Test_Y_ARMA02_Const(CheckArmaResultsMixin):
  265. @classmethod
  266. def setup_class(cls):
  267. endog = y_arma[:, 11]
  268. cls.res1 = ARMA(endog, order=(0, 2)).fit(trend="c", disp=-1)
  269. cls.res2 = results_arma.Y_arma02c()
  270. # cov_params and tvalues are off still but not as much vs. R
  271. class Test_Y_ARMA11_NoConst_CSS(CheckArmaResultsMixin):
  272. @classmethod
  273. def setup_class(cls):
  274. endog = y_arma[:, 0]
  275. cls.res1 = ARMA(endog, order=(1, 1)).fit(method="css", trend='nc',
  276. disp=-1)
  277. cls.res2 = results_arma.Y_arma11("css")
  278. cls.decimal_t = DECIMAL_1
  279. # better vs. R
  280. class Test_Y_ARMA14_NoConst_CSS(CheckArmaResultsMixin):
  281. @classmethod
  282. def setup_class(cls):
  283. endog = y_arma[:, 1]
  284. cls.res1 = ARMA(endog, order=(1, 4)).fit(method="css", trend='nc',
  285. disp=-1)
  286. cls.res2 = results_arma.Y_arma14("css")
  287. cls.decimal_fittedvalues = DECIMAL_3
  288. cls.decimal_resid = DECIMAL_3
  289. cls.decimal_t = DECIMAL_1
  290. # bse, etc. better vs. R
  291. # maroot is off because maparams is off a bit (adjust tolerance?)
  292. class Test_Y_ARMA41_NoConst_CSS(CheckArmaResultsMixin):
  293. @classmethod
  294. def setup_class(cls):
  295. endog = y_arma[:, 2]
  296. cls.res1 = ARMA(endog, order=(4, 1)).fit(method="css", trend='nc',
  297. disp=-1)
  298. cls.res2 = results_arma.Y_arma41("css")
  299. cls.decimal_t = DECIMAL_1
  300. cls.decimal_pvalues = 0
  301. cls.decimal_cov_params = DECIMAL_3
  302. cls.decimal_maroots = DECIMAL_1
  303. # same notes as above
  304. class Test_Y_ARMA22_NoConst_CSS(CheckArmaResultsMixin):
  305. @classmethod
  306. def setup_class(cls):
  307. endog = y_arma[:, 3]
  308. cls.res1 = ARMA(endog, order=(2, 2)).fit(method="css", trend='nc',
  309. disp=-1)
  310. cls.res2 = results_arma.Y_arma22("css")
  311. cls.decimal_t = DECIMAL_1
  312. cls.decimal_resid = DECIMAL_3
  313. cls.decimal_pvalues = DECIMAL_1
  314. cls.decimal_fittedvalues = DECIMAL_3
  315. # NOTE: gretl just uses least squares for AR CSS
  316. # so BIC, etc. is
  317. # -2*res1.llf + np.log(nobs)*(res1.q+res1.p+res1.k)
  318. # with no adjustment for p and no extra sigma estimate
  319. # NOTE: so our tests use x-12 arima results which agree with us and are
  320. # consistent with the rest of the models
  321. class Test_Y_ARMA50_NoConst_CSS(CheckArmaResultsMixin):
  322. @classmethod
  323. def setup_class(cls):
  324. endog = y_arma[:, 4]
  325. cls.res1 = ARMA(endog, order=(5, 0)).fit(method="css", trend='nc',
  326. disp=-1)
  327. cls.res2 = results_arma.Y_arma50("css")
  328. cls.decimal_t = 0
  329. cls.decimal_llf = DECIMAL_1 # looks like rounding error?
  330. class Test_Y_ARMA02_NoConst_CSS(CheckArmaResultsMixin):
  331. @classmethod
  332. def setup_class(cls):
  333. endog = y_arma[:, 5]
  334. cls.res1 = ARMA(endog, order=(0, 2)).fit(method="css", trend='nc',
  335. disp=-1)
  336. cls.res2 = results_arma.Y_arma02("css")
  337. # NOTE: our results are close to --x-12-arima option and R
  338. class Test_Y_ARMA11_Const_CSS(CheckArmaResultsMixin):
  339. @classmethod
  340. def setup_class(cls):
  341. endog = y_arma[:, 6]
  342. cls.res1 = ARMA(endog, order=(1, 1)).fit(trend="c", method="css",
  343. disp=-1)
  344. cls.res2 = results_arma.Y_arma11c("css")
  345. cls.decimal_params = DECIMAL_3
  346. cls.decimal_cov_params = DECIMAL_3
  347. cls.decimal_t = DECIMAL_1
  348. class Test_Y_ARMA14_Const_CSS(CheckArmaResultsMixin):
  349. @classmethod
  350. def setup_class(cls):
  351. endog = y_arma[:, 7]
  352. cls.res1 = ARMA(endog, order=(1, 4)).fit(trend="c", method="css",
  353. disp=-1)
  354. cls.res2 = results_arma.Y_arma14c("css")
  355. cls.decimal_t = DECIMAL_1
  356. cls.decimal_pvalues = DECIMAL_1
  357. class Test_Y_ARMA41_Const_CSS(CheckArmaResultsMixin):
  358. @classmethod
  359. def setup_class(cls):
  360. endog = y_arma[:, 8]
  361. cls.res1 = ARMA(endog, order=(4, 1)).fit(trend="c", method="css",
  362. disp=-1)
  363. cls.res2 = results_arma.Y_arma41c("css")
  364. cls.decimal_t = DECIMAL_1
  365. cls.decimal_cov_params = DECIMAL_1
  366. cls.decimal_maroots = DECIMAL_3
  367. cls.decimal_bse = DECIMAL_1
  368. class Test_Y_ARMA22_Const_CSS(CheckArmaResultsMixin):
  369. @classmethod
  370. def setup_class(cls):
  371. endog = y_arma[:, 9]
  372. cls.res1 = ARMA(endog, order=(2, 2)).fit(trend="c", method="css",
  373. disp=-1)
  374. cls.res2 = results_arma.Y_arma22c("css")
  375. cls.decimal_t = 0
  376. cls.decimal_pvalues = DECIMAL_1
  377. class Test_Y_ARMA50_Const_CSS(CheckArmaResultsMixin):
  378. @classmethod
  379. def setup_class(cls):
  380. endog = y_arma[:, 10]
  381. cls.res1 = ARMA(endog, order=(5, 0)).fit(trend="c", method="css",
  382. disp=-1)
  383. cls.res2 = results_arma.Y_arma50c("css")
  384. cls.decimal_t = DECIMAL_1
  385. cls.decimal_params = DECIMAL_3
  386. cls.decimal_cov_params = DECIMAL_2
  387. class Test_Y_ARMA02_Const_CSS(CheckArmaResultsMixin):
  388. @classmethod
  389. def setup_class(cls):
  390. endog = y_arma[:, 11]
  391. cls.res1 = ARMA(endog, order=(0, 2)).fit(trend="c", method="css",
  392. disp=-1)
  393. cls.res2 = results_arma.Y_arma02c("css")
  394. def test_reset_trend_error():
  395. endog = y_arma[:, 0]
  396. mod = ARMA(endog, order=(1, 1))
  397. mod.fit(trend="c", disp=-1)
  398. with pytest.raises(RuntimeError):
  399. mod.fit(trend="nc", disp=-1)
  400. @pytest.mark.slow
  401. def test_start_params_bug():
  402. data = np.array([1368., 1187, 1090, 1439, 2362, 2783, 2869, 2512, 1804,
  403. 1544, 1028, 869, 1737, 2055, 1947, 1618, 1196, 867, 997,
  404. 1862, 2525,
  405. 3250, 4023, 4018, 3585, 3004, 2500, 2441, 2749, 2466,
  406. 2157, 1847, 1463,
  407. 1146, 851, 993, 1448, 1719, 1709, 1455, 1950, 1763, 2075,
  408. 2343, 3570,
  409. 4690, 3700, 2339, 1679, 1466, 998, 853, 835, 922, 851,
  410. 1125, 1299, 1105,
  411. 860, 701, 689, 774, 582, 419, 846, 1132, 902, 1058, 1341,
  412. 1551, 1167,
  413. 975, 786, 759, 751, 649, 876, 720, 498, 553, 459, 543,
  414. 447, 415, 377,
  415. 373, 324, 320, 306, 259, 220, 342, 558, 825, 994, 1267,
  416. 1473, 1601,
  417. 1896, 1890, 2012, 2198, 2393, 2825, 3411, 3406, 2464,
  418. 2891, 3685, 3638,
  419. 3746, 3373, 3190, 2681, 2846, 4129, 5054, 5002, 4801,
  420. 4934, 4903, 4713,
  421. 4745, 4736, 4622, 4642, 4478, 4510, 4758, 4457, 4356,
  422. 4170, 4658, 4546,
  423. 4402, 4183, 3574, 2586, 3326, 3948, 3983, 3997, 4422,
  424. 4496, 4276, 3467,
  425. 2753, 2582, 2921, 2768, 2789, 2824, 2482, 2773, 3005,
  426. 3641, 3699, 3774,
  427. 3698, 3628, 3180, 3306, 2841, 2014, 1910, 2560, 2980,
  428. 3012, 3210, 3457,
  429. 3158, 3344, 3609, 3327, 2913, 2264, 2326, 2596, 2225,
  430. 1767, 1190, 792,
  431. 669, 589, 496, 354, 246, 250, 323, 495, 924, 1536, 2081,
  432. 2660, 2814, 2992,
  433. 3115, 2962, 2272, 2151, 1889, 1481, 955, 631, 288, 103,
  434. 60, 82, 107, 185,
  435. 618, 1526, 2046, 2348, 2584, 2600, 2515, 2345, 2351, 2355,
  436. 2409, 2449,
  437. 2645, 2918, 3187, 2888, 2610, 2740, 2526, 2383, 2936,
  438. 2968, 2635, 2617,
  439. 2790, 3906, 4018, 4797, 4919, 4942, 4656, 4444, 3898,
  440. 3908, 3678, 3605,
  441. 3186, 2139, 2002, 1559, 1235, 1183, 1096, 673, 389, 223,
  442. 352, 308, 365,
  443. 525, 779, 894, 901, 1025, 1047, 981, 902, 759, 569, 519,
  444. 408, 263, 156,
  445. 72, 49, 31, 41, 192, 423, 492, 552, 564, 723, 921, 1525,
  446. 2768, 3531, 3824,
  447. 3835, 4294, 4533, 4173, 4221, 4064, 4641, 4685, 4026,
  448. 4323, 4585, 4836,
  449. 4822, 4631, 4614, 4326, 4790, 4736, 4104, 5099, 5154,
  450. 5121, 5384, 5274,
  451. 5225, 4899, 5382, 5295, 5349, 4977, 4597, 4069, 3733,
  452. 3439, 3052, 2626,
  453. 1939, 1064, 713, 916, 832, 658, 817, 921, 772, 764, 824,
  454. 967, 1127, 1153,
  455. 824, 912, 957, 990, 1218, 1684, 2030, 2119, 2233, 2657,
  456. 2652, 2682, 2498,
  457. 2429, 2346, 2298, 2129, 1829, 1816, 1225, 1010, 748, 627,
  458. 469, 576, 532,
  459. 475, 582, 641, 605, 699, 680, 714, 670, 666, 636, 672,
  460. 679, 446, 248, 134,
  461. 160, 178, 286, 413, 676, 1025, 1159, 952, 1398, 1833,
  462. 2045, 2072, 1798,
  463. 1799, 1358, 727, 353, 347, 844, 1377, 1829, 2118, 2272,
  464. 2745, 4263, 4314,
  465. 4530, 4354, 4645, 4547, 5391, 4855, 4739, 4520, 4573,
  466. 4305, 4196, 3773,
  467. 3368, 2596, 2596, 2305, 2756, 3747, 4078, 3415, 2369,
  468. 2210, 2316, 2263,
  469. 2672, 3571, 4131, 4167, 4077, 3924, 3738, 3712, 3510,
  470. 3182, 3179, 2951,
  471. 2453, 2078, 1999, 2486, 2581, 1891, 1997, 1366, 1294,
  472. 1536, 2794, 3211,
  473. 3242, 3406, 3121, 2425, 2016, 1787, 1508, 1304, 1060,
  474. 1342, 1589, 2361,
  475. 3452, 2659, 2857, 3255, 3322, 2852, 2964, 3132, 3033,
  476. 2931, 2636, 2818,
  477. 3310, 3396, 3179, 3232, 3543, 3759, 3503, 3758, 3658,
  478. 3425, 3053, 2620,
  479. 1837, 923, 712, 1054, 1376, 1556, 1498, 1523, 1088, 728,
  480. 890, 1413, 2524,
  481. 3295, 4097, 3993, 4116, 3874, 4074, 4142, 3975, 3908,
  482. 3907, 3918, 3755,
  483. 3648, 3778, 4293, 4385, 4360, 4352, 4528, 4365, 3846,
  484. 4098, 3860, 3230,
  485. 2820, 2916, 3201, 3721, 3397, 3055, 2141, 1623, 1825,
  486. 1716, 2232, 2939,
  487. 3735, 4838, 4560, 4307, 4975, 5173, 4859, 5268, 4992,
  488. 5100, 5070, 5270,
  489. 4760, 5135, 5059, 4682, 4492, 4933, 4737, 4611, 4634,
  490. 4789, 4811, 4379,
  491. 4689, 4284, 4191, 3313, 2770, 2543, 3105, 2967, 2420,
  492. 1996, 2247, 2564,
  493. 2726, 3021, 3427, 3509, 3759, 3324, 2988, 2849, 2340,
  494. 2443, 2364, 1252,
  495. 623, 742, 867, 684, 488, 348, 241, 187, 279, 355, 423,
  496. 678, 1375, 1497,
  497. 1434, 2116, 2411, 1929, 1628, 1635, 1609, 1757, 2090,
  498. 2085, 1790, 1846,
  499. 2038, 2360, 2342, 2401, 2920, 3030, 3132, 4385, 5483,
  500. 5865, 5595, 5485,
  501. 5727, 5553, 5560, 5233, 5478, 5159, 5155, 5312, 5079,
  502. 4510, 4628, 4535,
  503. 3656, 3698, 3443, 3146, 2562, 2304, 2181, 2293, 1950,
  504. 1930, 2197, 2796,
  505. 3441, 3649, 3815, 2850, 4005, 5305, 5550, 5641, 4717,
  506. 5131, 2831, 3518,
  507. 3354, 3115, 3515, 3552, 3244, 3658, 4407, 4935, 4299,
  508. 3166, 3335, 2728,
  509. 2488, 2573, 2002, 1717, 1645, 1977, 2049, 2125, 2376,
  510. 2551, 2578, 2629,
  511. 2750, 3150, 3699, 4062, 3959, 3264, 2671, 2205, 2128,
  512. 2133, 2095, 1964,
  513. 2006, 2074, 2201, 2506, 2449, 2465, 2064, 1446, 1382, 983,
  514. 898, 489, 319,
  515. 383, 332, 276, 224, 144, 101, 232, 429, 597, 750, 908,
  516. 960, 1076, 951,
  517. 1062, 1183, 1404, 1391, 1419, 1497, 1267, 963, 682, 777,
  518. 906, 1149, 1439,
  519. 1600, 1876, 1885, 1962, 2280, 2711, 2591, 2411])
  520. with warnings.catch_warnings():
  521. warnings.simplefilter("ignore")
  522. ARMA(data, order=(4, 1)).fit(start_ar_lags=5, disp=-1)
  523. class Test_ARIMA101(CheckArmaResultsMixin):
  524. @classmethod
  525. def setup_class(cls):
  526. endog = y_arma[:, 6]
  527. cls.res1 = ARIMA(endog, (1, 0, 1)).fit(trend="c", disp=-1)
  528. (cls.res1.forecast_res, cls.res1.forecast_err,
  529. confint) = cls.res1.forecast(10)
  530. cls.res2 = results_arma.Y_arma11c()
  531. cls.res2.k_diff = 0
  532. cls.res2.k_ar = 1
  533. cls.res2.k_ma = 1
  534. class Test_ARIMA111(CheckArimaResultsMixin, CheckForecastMixin,
  535. CheckDynamicForecastMixin):
  536. @classmethod
  537. def setup_class(cls):
  538. cpi = load_macrodata_pandas().data['cpi'].values
  539. cls.res1 = ARIMA(cpi, (1, 1, 1)).fit(disp=-1)
  540. cls.res2 = results_arima.ARIMA111()
  541. # make sure endog names changes to D.cpi
  542. cls.decimal_llf = 3
  543. cls.decimal_aic = 3
  544. cls.decimal_bic = 3
  545. # TODO: why has dec_cov_params changed, used to be better
  546. cls.decimal_cov_params = 2
  547. cls.decimal_t = 0
  548. (cls.res1.forecast_res,
  549. cls.res1.forecast_err,
  550. conf_int) = cls.res1.forecast(25)
  551. # TODO: fix the indexing for the end here, I do not think this is right
  552. # if we're going to treat it like indexing
  553. # the forecast from 2005Q1 through 2009Q4 is indices
  554. # 184 through 227 not 226
  555. # note that the first one counts in the count so 164 + 64 is 65
  556. # predictions
  557. cls.res1.forecast_res_dyn = cls.res1.predict(start=164, end=164 + 63,
  558. typ='levels',
  559. dynamic=True)
  560. def test_freq(self):
  561. assert_almost_equal(self.res1.arfreq, [0.0000], 4)
  562. assert_almost_equal(self.res1.mafreq, [0.0000], 4)
  563. class Test_ARIMA111CSS(CheckArimaResultsMixin, CheckForecastMixin,
  564. CheckDynamicForecastMixin):
  565. @classmethod
  566. def setup_class(cls):
  567. cpi = load_macrodata_pandas().data['cpi'].values
  568. cls.res1 = ARIMA(cpi, (1, 1, 1)).fit(disp=-1, method='css')
  569. cls.res2 = results_arima.ARIMA111(method='css')
  570. cls.res2.fittedvalues = - cpi[1:-1] + cls.res2.linear
  571. # make sure endog names changes to D.cpi
  572. (cls.res1.forecast_res,
  573. cls.res1.forecast_err,
  574. conf_int) = cls.res1.forecast(25)
  575. cls.decimal_forecast = 2
  576. cls.decimal_forecast_dyn = 2
  577. cls.decimal_forecasterr = 3
  578. cls.res1.forecast_res_dyn = cls.res1.predict(start=164, end=164 + 63,
  579. typ='levels',
  580. dynamic=True)
  581. # precisions
  582. cls.decimal_arroots = 3
  583. cls.decimal_cov_params = 3
  584. cls.decimal_hqic = 3
  585. cls.decimal_maroots = 3
  586. cls.decimal_t = 1
  587. cls.decimal_fittedvalues = 2 # because of rounding when copying
  588. cls.decimal_resid = 2
  589. cls.decimal_predict_levels = DECIMAL_2
  590. class Test_ARIMA112CSS(CheckArimaResultsMixin):
  591. @classmethod
  592. def setup_class(cls):
  593. cpi = load_macrodata_pandas().data['cpi'].values
  594. cls.res1 = ARIMA(cpi, (1, 1, 2)).fit(disp=-1, method='css',
  595. start_params=[.905322, -.692425,
  596. 1.07366,
  597. 0.172024])
  598. cls.res2 = results_arima.ARIMA112(method='css')
  599. cls.res2.fittedvalues = - cpi[1:-1] + cls.res2.linear
  600. # make sure endog names changes to D.cpi
  601. cls.decimal_llf = 3
  602. cls.decimal_aic = 3
  603. cls.decimal_bic = 3
  604. # TODO: fix the indexing for the end here, I do not think this is right
  605. # if we're going to treat it like indexing
  606. # the forecast from 2005Q1 through 2009Q4 is indices
  607. # 184 through 227 not 226
  608. # note that the first one counts in the count so 164 + 64 is 65
  609. # predictions
  610. # cls.res1.forecast_res_dyn = self.predict(start=164, end=164+63,
  611. # typ='levels', dynamic=True)
  612. # since we got from gretl do not have linear prediction in differences
  613. cls.decimal_arroots = 3
  614. cls.decimal_maroots = 2
  615. cls.decimal_t = 1
  616. cls.decimal_resid = 2
  617. cls.decimal_fittedvalues = 3
  618. cls.decimal_predict_levels = DECIMAL_3
  619. def test_freq(self):
  620. assert_almost_equal(self.res1.arfreq, [0.5000], 4)
  621. assert_almost_equal(self.res1.mafreq, [0.5000, 0.5000], 4)
  622. def test_arima_predict_mle_dates():
  623. cpi = load_macrodata_pandas().data['cpi'].values
  624. res1 = ARIMA(cpi, (4, 1, 1), dates=cpi_dates, freq='Q').fit(disp=-1)
  625. file_path = os.path.join(current_path, 'results',
  626. 'results_arima_forecasts_all_mle.csv')
  627. with open(file_path, "rb") as test_data:
  628. arima_forecasts = np.genfromtxt(test_data, delimiter=",",
  629. skip_header=1, dtype=float)
  630. fc = arima_forecasts[:, 0]
  631. fcdyn = arima_forecasts[:, 1]
  632. fcdyn2 = arima_forecasts[:, 2]
  633. start, end = 2, 51
  634. fv = res1.predict('1959Q3', '1971Q4', typ='levels')
  635. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  636. assert_equal(res1.data.predict_dates, cpi_dates[start:end + 1])
  637. start, end = 202, 227
  638. fv = res1.predict('2009Q3', '2015Q4', typ='levels')
  639. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  640. assert_equal(res1.data.predict_dates, cpi_predict_dates)
  641. # make sure dynamic works
  642. start, end = '1960q2', '1971q4'
  643. fv = res1.predict(start, end, dynamic=True, typ='levels')
  644. assert_almost_equal(fv, fcdyn[5:51 + 1], DECIMAL_4)
  645. start, end = '1965q1', '2015q4'
  646. fv = res1.predict(start, end, dynamic=True, typ='levels')
  647. assert_almost_equal(fv, fcdyn2[24:227 + 1], DECIMAL_4)
  648. def test_arma_predict_mle_dates():
  649. from statsmodels.datasets.sunspots import load_pandas
  650. sunspots = load_pandas().data['SUNACTIVITY'].values
  651. mod = ARMA(sunspots, (9, 0), dates=sun_dates, freq='A')
  652. mod.method = 'mle'
  653. assert_raises(ValueError, mod._get_prediction_index, '1701', '1751', True)
  654. start, end = 2, 51
  655. mod._get_prediction_index('1702', '1751', False)
  656. assert_equal(mod.data.predict_dates, sun_dates[start:end + 1])
  657. start, end = 308, 333
  658. mod._get_prediction_index('2008', '2033', False)
  659. assert_equal(mod.data.predict_dates, sun_predict_dates)
  660. def test_arima_predict_css_dates():
  661. cpi = load_macrodata_pandas().data['cpi'].values
  662. res1 = ARIMA(cpi, (4, 1, 1), dates=cpi_dates, freq='Q').fit(disp=-1,
  663. method='css',
  664. trend='nc')
  665. params = np.array([1.231272508473910,
  666. -0.282516097759915,
  667. 0.170052755782440,
  668. -0.118203728504945,
  669. -0.938783134717947])
  670. file_path = os.path.join(current_path, 'results',
  671. 'results_arima_forecasts_all_css.csv')
  672. with open(file_path, "rb") as test_data:
  673. arima_forecasts = np.genfromtxt(test_data, delimiter=",",
  674. skip_header=1, dtype=float)
  675. fc = arima_forecasts[:, 0]
  676. fcdyn = arima_forecasts[:, 1]
  677. fcdyn2 = arima_forecasts[:, 2]
  678. start, end = 5, 51
  679. fv = res1.model.predict(params, '1960Q2', '1971Q4', typ='levels')
  680. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  681. assert_equal(res1.data.predict_dates, cpi_dates[start:end + 1])
  682. start, end = 202, 227
  683. fv = res1.model.predict(params, '2009Q3', '2015Q4', typ='levels')
  684. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  685. assert_equal(res1.data.predict_dates, cpi_predict_dates)
  686. # make sure dynamic works
  687. start, end = 5, 51
  688. fv = res1.model.predict(params, '1960Q2', '1971Q4', typ='levels',
  689. dynamic=True)
  690. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  691. start, end = '1965q1', '2015q4'
  692. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  693. assert_almost_equal(fv, fcdyn2[24:227 + 1], DECIMAL_4)
  694. def test_arma_predict_css_dates():
  695. from statsmodels.datasets.sunspots import load_pandas
  696. sunspots = load_pandas().data['SUNACTIVITY'].values
  697. mod = ARMA(sunspots, (9, 0), dates=sun_dates, freq='A')
  698. mod.method = 'css'
  699. assert_raises(ValueError, mod._get_prediction_index, '1701', '1751', False)
  700. def test_arima_predict_mle():
  701. cpi = load_macrodata_pandas().data['cpi'].values
  702. res1 = ARIMA(cpi, (4, 1, 1)).fit(disp=-1)
  703. # fit the model so that we get correct endog length but use
  704. file_path = os.path.join(current_path, 'results',
  705. 'results_arima_forecasts_all_mle.csv')
  706. with open(file_path, "rb") as test_data:
  707. arima_forecasts = np.genfromtxt(test_data, delimiter=",",
  708. skip_header=1, dtype=float)
  709. fc = arima_forecasts[:, 0]
  710. fcdyn = arima_forecasts[:, 1]
  711. fcdyn2 = arima_forecasts[:, 2]
  712. fcdyn3 = arima_forecasts[:, 3]
  713. fcdyn4 = arima_forecasts[:, 4]
  714. # 0 indicates the first sample-observation below
  715. # ie., the index after the pre-sample, these are also differenced once
  716. # so the indices are moved back once from the cpi in levels
  717. # start < p, end <p 1959q2 - 1959q4
  718. start, end = 1, 3
  719. fv = res1.predict(start, end, typ='levels')
  720. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  721. # start < p, end 0 1959q3 - 1960q1
  722. start, end = 2, 4
  723. fv = res1.predict(start, end, typ='levels')
  724. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  725. # start < p, end >0 1959q3 - 1971q4
  726. start, end = 2, 51
  727. fv = res1.predict(start, end, typ='levels')
  728. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  729. # start < p, end nobs 1959q3 - 2009q3
  730. start, end = 2, 202
  731. fv = res1.predict(start, end, typ='levels')
  732. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  733. # start < p, end >nobs 1959q3 - 2015q4
  734. start, end = 2, 227
  735. fv = res1.predict(start, end, typ='levels')
  736. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  737. # start 0, end >0 1960q1 - 1971q4
  738. start, end = 4, 51
  739. fv = res1.predict(start, end, typ='levels')
  740. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  741. # start 0, end nobs 1960q1 - 2009q3
  742. start, end = 4, 202
  743. fv = res1.predict(start, end, typ='levels')
  744. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  745. # start 0, end >nobs 1960q1 - 2015q4
  746. start, end = 4, 227
  747. fv = res1.predict(start, end, typ='levels')
  748. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  749. # start >p, end >0 1965q1 - 1971q4
  750. start, end = 24, 51
  751. fv = res1.predict(start, end, typ='levels')
  752. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  753. # start >p, end nobs 1965q1 - 2009q3
  754. start, end = 24, 202
  755. fv = res1.predict(start, end, typ='levels')
  756. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  757. # start >p, end >nobs 1965q1 - 2015q4
  758. start, end = 24, 227
  759. fv = res1.predict(start, end, typ='levels')
  760. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  761. # start nobs, end nobs 2009q3 - 2009q3
  762. start, end = 202, 202
  763. fv = res1.predict(start, end, typ='levels')
  764. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_3)
  765. # start nobs, end >nobs 2009q3 - 2015q4
  766. start, end = 202, 227
  767. fv = res1.predict(start, end, typ='levels')
  768. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_3)
  769. # start >nobs, end >nobs 2009q4 - 2015q4
  770. start, end = 203, 227
  771. fv = res1.predict(start, end, typ='levels')
  772. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  773. # defaults
  774. start, end = None, None
  775. fv = res1.predict(start, end, typ='levels')
  776. assert_almost_equal(fv, fc[1:203], DECIMAL_4)
  777. # Dynamic
  778. # start < p, end <p 1959q2 - 1959q4
  779. start, end = 1, 3
  780. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  781. fv = res1.predict(start, end, dynamic=True, typ='levels')
  782. # start < p, end 0 1959q3 - 1960q1
  783. start, end = 2, 4
  784. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  785. res1.predict(start, end, dynamic=True, typ='levels')
  786. # start < p, end >0 1959q3 - 1971q4
  787. start, end = 2, 51
  788. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  789. res1.predict(start, end, dynamic=True, typ='levels')
  790. # start < p, end nobs 1959q3 - 2009q3
  791. start, end = 2, 202
  792. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  793. res1.predict(start, end, dynamic=True, typ='levels')
  794. # start < p, end >nobs 1959q3 - 2015q4
  795. start, end = 2, 227
  796. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  797. res1.predict(start, end, dynamic=True, typ='levels')
  798. # start 0, end >0 1960q1 - 1971q4
  799. start, end = 5, 51
  800. fv = res1.predict(start, end, dynamic=True, typ='levels')
  801. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  802. # start 0, end nobs 1960q1 - 2009q3
  803. start, end = 5, 202
  804. fv = res1.predict(start, end, dynamic=True, typ='levels')
  805. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  806. # start 0, end >nobs 1960q1 - 2015q4
  807. start, end = 5, 227
  808. fv = res1.predict(start, end, dynamic=True, typ='levels')
  809. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  810. # start >p, end >0 1965q1 - 1971q4
  811. start, end = 24, 51
  812. fv = res1.predict(start, end, dynamic=True, typ='levels')
  813. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  814. # start >p, end nobs 1965q1 - 2009q3
  815. start, end = 24, 202
  816. fv = res1.predict(start, end, dynamic=True, typ='levels')
  817. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  818. # start >p, end >nobs 1965q1 - 2015q4
  819. start, end = 24, 227
  820. fv = res1.predict(start, end, dynamic=True, typ='levels')
  821. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  822. # start nobs, end nobs 2009q3 - 2009q3
  823. start, end = 202, 202
  824. fv = res1.predict(start, end, dynamic=True, typ='levels')
  825. assert_almost_equal(fv, fcdyn3[start:end + 1], DECIMAL_4)
  826. # start nobs, end >nobs 2009q3 - 2015q4
  827. start, end = 202, 227
  828. fv = res1.predict(start, end, dynamic=True, typ='levels')
  829. assert_almost_equal(fv, fcdyn3[start:end + 1], DECIMAL_4)
  830. # start >nobs, end >nobs 2009q4 - 2015q4
  831. start, end = 203, 227
  832. fv = res1.predict(start, end, dynamic=True, typ='levels')
  833. assert_almost_equal(fv, fcdyn4[start:end + 1], DECIMAL_4)
  834. # defaults
  835. start, end = None, None
  836. fv = res1.predict(start, end, dynamic=True, typ='levels')
  837. assert_almost_equal(fv, fcdyn[5:203], DECIMAL_4)
  838. def _check_start(model, given, expected, dynamic):
  839. start, _, _, _ = model._get_prediction_index(given, None, dynamic)
  840. assert_equal(start, expected)
  841. def _check_end(model, given, end_expect, out_of_sample_expect):
  842. _, end, out_of_sample, _ = model._get_prediction_index(None, given, False)
  843. assert_equal((end, out_of_sample), (end_expect, out_of_sample_expect))
  844. def test_arma_predict_indices():
  845. from statsmodels.datasets.sunspots import load_pandas
  846. sunspots = load_pandas().data['SUNACTIVITY'].values
  847. model = ARMA(sunspots, (9, 0), dates=sun_dates, freq='A')
  848. model.method = 'mle'
  849. # raises - pre-sample + dynamic
  850. with pytest.raises(ValueError):
  851. model._get_prediction_index(0, None, True)
  852. with pytest.raises(ValueError):
  853. model._get_prediction_index(8, None, True)
  854. with pytest.raises(ValueError):
  855. model._get_prediction_index('1700', None, True)
  856. with pytest.raises(ValueError):
  857. model._get_prediction_index('1708', None, True)
  858. # raises - start out of sample
  859. # works - in-sample
  860. # None
  861. start_test_cases = [
  862. # given, expected, dynamic
  863. (None, 9, True),
  864. # all start get moved back by k_diff
  865. (9, 9, True),
  866. (10, 10, True),
  867. # what about end of sample start - last value is first
  868. # forecast
  869. (309, 309, True),
  870. (308, 308, True),
  871. (0, 0, False),
  872. (1, 1, False),
  873. (4, 4, False),
  874. # all start get moved back by k_diff
  875. ('1709', 9, True),
  876. ('1710', 10, True),
  877. # what about end of sample start - last value is first
  878. # forecast
  879. ('2008', 308, True),
  880. ('2009', 309, True),
  881. ('1700', 0, False),
  882. ('1708', 8, False),
  883. ('1709', 9, False),
  884. ]
  885. for case in start_test_cases:
  886. _check_start(*((model,) + case))
  887. # the length of sunspot is 309, so last index is 208
  888. end_test_cases = [(None, 308, 0),
  889. (307, 307, 0),
  890. (308, 308, 0),
  891. (309, 308, 1),
  892. (312, 308, 4),
  893. (51, 51, 0),
  894. (333, 308, 25),
  895. ('2007', 307, 0),
  896. ('2008', 308, 0),
  897. ('2009', 308, 1),
  898. ('2012', 308, 4),
  899. ('1815', 115, 0),
  900. ('2033', 308, 25),
  901. ]
  902. for case in end_test_cases:
  903. _check_end(*((model,) + case))
  904. def test_arima_predict_indices():
  905. cpi = load_macrodata_pandas().data['cpi'].values
  906. model = ARIMA(cpi, (4, 1, 1), dates=cpi_dates, freq='Q')
  907. model.method = 'mle'
  908. # starting indices
  909. # raises - pre-sample + dynamic
  910. with pytest.raises(ValueError):
  911. model._get_prediction_index(0, None, True)
  912. with pytest.raises(ValueError):
  913. model._get_prediction_index(4, None, True)
  914. with pytest.raises(KeyError):
  915. model._get_prediction_index('1959Q1', None, True)
  916. with pytest.raises(ValueError):
  917. model._get_prediction_index('1960Q1', None, True)
  918. # raises - index differenced away
  919. with pytest.raises(ValueError):
  920. model._get_prediction_index(0, None, False)
  921. with pytest.raises(KeyError):
  922. model._get_prediction_index('1959Q1', None, False)
  923. # raises - start out of sample
  924. # works - in-sample
  925. # None
  926. start_test_cases = [
  927. # given, expected, dynamic
  928. (None, 4, True),
  929. # all start get moved back by k_diff
  930. (5, 4, True),
  931. (6, 5, True),
  932. # what about end of sample start - last value is first
  933. # forecast
  934. (203, 202, True),
  935. (1, 0, False),
  936. (4, 3, False),
  937. (5, 4, False),
  938. # all start get moved back by k_diff
  939. ('1960Q2', 4, True),
  940. ('1960Q3', 5, True),
  941. # what about end of sample start - last value is first
  942. # forecast
  943. ('2009Q4', 202, True),
  944. ('1959Q2', 0, False),
  945. ('1960Q1', 3, False),
  946. ('1960Q2', 4, False),
  947. ]
  948. for case in start_test_cases:
  949. _check_start(*((model,) + case))
  950. # TODO: make sure dates are passing through unmolested
  951. # the length of diff(cpi) is 202, so last index is 201
  952. end_test_cases = [(None, 201, 0),
  953. (201, 200, 0),
  954. (202, 201, 0),
  955. (203, 201, 1),
  956. (204, 201, 2),
  957. (51, 50, 0),
  958. (164 + 63, 201, 25),
  959. ('2009Q2', 200, 0),
  960. ('2009Q3', 201, 0),
  961. ('2009Q4', 201, 1),
  962. ('2010Q1', 201, 2),
  963. ('1971Q4', 50, 0),
  964. ('2015Q4', 201, 25),
  965. ]
  966. for case in end_test_cases:
  967. _check_end(*((model,) + case))
  968. # check higher k_diff
  969. # model.k_diff = 2
  970. model = ARIMA(cpi, (4, 2, 1), dates=cpi_dates, freq='Q')
  971. model.method = 'mle'
  972. # raises - pre-sample + dynamic
  973. assert_raises(ValueError, model._get_prediction_index, 0, None, True)
  974. assert_raises(ValueError, model._get_prediction_index, 5, None, True)
  975. assert_raises(KeyError, model._get_prediction_index,
  976. '1959Q1', None, True)
  977. assert_raises(ValueError, model._get_prediction_index,
  978. '1960Q1', None, True)
  979. # raises - index differenced away
  980. assert_raises(ValueError, model._get_prediction_index, 1, None, False)
  981. assert_raises(KeyError, model._get_prediction_index,
  982. '1959Q2', None, False)
  983. start_test_cases = [(None, 4, True),
  984. # all start get moved back by k_diff
  985. (6, 4, True),
  986. # what about end of sample start - last value is first
  987. # forecast
  988. (203, 201, True),
  989. (2, 0, False),
  990. (4, 2, False),
  991. (5, 3, False),
  992. ('1960Q3', 4, True),
  993. # what about end of sample start - last value is first
  994. # forecast
  995. ('2009Q4', 201, True),
  996. ('2009Q4', 201, True),
  997. ('1959Q3', 0, False),
  998. ('1960Q1', 2, False),
  999. ('1960Q2', 3, False),
  1000. ]
  1001. for case in start_test_cases:
  1002. _check_start(*((model,) + case))
  1003. end_test_cases = [(None, 200, 0),
  1004. (201, 199, 0),
  1005. (202, 200, 0),
  1006. (203, 200, 1),
  1007. (204, 200, 2),
  1008. (51, 49, 0),
  1009. (164 + 63, 200, 25),
  1010. ('2009Q2', 199, 0),
  1011. ('2009Q3', 200, 0),
  1012. ('2009Q4', 200, 1),
  1013. ('2010Q1', 200, 2),
  1014. ('1971Q4', 49, 0),
  1015. ('2015Q4', 200, 25),
  1016. ]
  1017. for case in end_test_cases:
  1018. _check_end(*((model,) + case))
  1019. def test_arima_predict_indices_css():
  1020. cpi = load_macrodata_pandas().data['cpi'].values
  1021. # NOTE: Doing no-constant for now to kick the conditional exogenous
  1022. # issue 274 down the road
  1023. # go ahead and git the model to set up necessary variables
  1024. model = ARIMA(cpi, (4, 1, 1))
  1025. model.method = 'css'
  1026. assert_raises(ValueError, model._get_prediction_index, 0, None, False)
  1027. assert_raises(ValueError, model._get_prediction_index, 0, None, True)
  1028. assert_raises(ValueError, model._get_prediction_index, 2, None, False)
  1029. assert_raises(ValueError, model._get_prediction_index, 2, None, True)
  1030. def test_arima_predict_css():
  1031. cpi = load_macrodata_pandas().data['cpi'].values
  1032. # NOTE: Doing no-constant for now to kick the conditional exogenous
  1033. # issue 274 down the road
  1034. # go ahead and git the model to set up necessary variables
  1035. res1 = ARIMA(cpi, (4, 1, 1)).fit(disp=-1, method="css",
  1036. trend="nc")
  1037. # but use gretl parameters to predict to avoid precision problems
  1038. params = np.array([1.231272508473910,
  1039. -0.282516097759915,
  1040. 0.170052755782440,
  1041. -0.118203728504945,
  1042. -0.938783134717947])
  1043. file_path = os.path.join(current_path, 'results',
  1044. 'results_arima_forecasts_all_css.csv')
  1045. with open(file_path, "rb") as test_data:
  1046. arima_forecasts = np.genfromtxt(test_data, delimiter=",",
  1047. skip_header=1, dtype=float)
  1048. fc = arima_forecasts[:, 0]
  1049. fcdyn = arima_forecasts[:, 1]
  1050. fcdyn2 = arima_forecasts[:, 2]
  1051. fcdyn3 = arima_forecasts[:, 3]
  1052. fcdyn4 = arima_forecasts[:, 4]
  1053. start, end = 1, 3
  1054. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1055. res1.model.predict(params, start, end)
  1056. # start < p, end 0 1959q3 - 1960q1
  1057. start, end = 2, 4
  1058. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1059. res1.model.predict(params, start, end)
  1060. # start < p, end >0 1959q3 - 1971q4
  1061. start, end = 2, 51
  1062. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1063. res1.model.predict(params, start, end)
  1064. # start < p, end nobs 1959q3 - 2009q3
  1065. start, end = 2, 202
  1066. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1067. res1.model.predict(params, start, end)
  1068. # start < p, end >nobs 1959q3 - 2015q4
  1069. start, end = 2, 227
  1070. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1071. res1.model.predict(params, start, end)
  1072. # start 0, end >0 1960q1 - 1971q4
  1073. start, end = 5, 51
  1074. fv = res1.model.predict(params, start, end, typ='levels')
  1075. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1076. # start 0, end nobs 1960q1 - 2009q3
  1077. start, end = 5, 202
  1078. fv = res1.model.predict(params, start, end, typ='levels')
  1079. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1080. # start 0, end >nobs 1960q1 - 2015q4
  1081. # TODO: why detoriating precision?
  1082. fv = res1.model.predict(params, start, end, typ='levels')
  1083. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1084. # start >p, end >0 1965q1 - 1971q4
  1085. start, end = 24, 51
  1086. fv = res1.model.predict(params, start, end, typ='levels')
  1087. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1088. # start >p, end nobs 1965q1 - 2009q3
  1089. start, end = 24, 202
  1090. fv = res1.model.predict(params, start, end, typ='levels')
  1091. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1092. # start >p, end >nobs 1965q1 - 2015q4
  1093. start, end = 24, 227
  1094. fv = res1.model.predict(params, start, end, typ='levels')
  1095. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1096. # start nobs, end nobs 2009q3 - 2009q3
  1097. start, end = 202, 202
  1098. fv = res1.model.predict(params, start, end, typ='levels')
  1099. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1100. # start nobs, end >nobs 2009q3 - 2015q4
  1101. start, end = 202, 227
  1102. fv = res1.model.predict(params, start, end, typ='levels')
  1103. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1104. # start >nobs, end >nobs 2009q4 - 2015q4
  1105. start, end = 203, 227
  1106. fv = res1.model.predict(params, start, end, typ='levels')
  1107. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1108. # defaults
  1109. start, end = None, None
  1110. fv = res1.model.predict(params, start, end, typ='levels')
  1111. assert_almost_equal(fv, fc[5:203], DECIMAL_4)
  1112. # Dynamic
  1113. # start < p, end <p 1959q2 - 1959q4
  1114. start, end = 1, 3
  1115. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1116. res1.predict(start, end, dynamic=True)
  1117. # start < p, end 0 1959q3 - 1960q1
  1118. start, end = 2, 4
  1119. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1120. res1.predict(start, end, dynamic=True)
  1121. # start < p, end >0 1959q3 - 1971q4
  1122. start, end = 2, 51
  1123. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1124. res1.predict(start, end, dynamic=True)
  1125. # start < p, end nobs 1959q3 - 2009q3
  1126. start, end = 2, 202
  1127. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1128. res1.predict(start, end, dynamic=True)
  1129. # start < p, end >nobs 1959q3 - 2015q4
  1130. start, end = 2, 227
  1131. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1132. res1.predict(start, end, dynamic=True)
  1133. # start 0, end >0 1960q1 - 1971q4
  1134. start, end = 5, 51
  1135. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1136. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1137. # start 0, end nobs 1960q1 - 2009q3
  1138. start, end = 5, 202
  1139. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1140. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1141. # start 0, end >nobs 1960q1 - 2015q4
  1142. start, end = 5, 227
  1143. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1144. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1145. # start >p, end >0 1965q1 - 1971q4
  1146. start, end = 24, 51
  1147. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1148. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1149. # start >p, end nobs 1965q1 - 2009q3
  1150. start, end = 24, 202
  1151. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1152. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1153. # start >p, end >nobs 1965q1 - 2015q4
  1154. start, end = 24, 227
  1155. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1156. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1157. # start nobs, end nobs 2009q3 - 2009q3
  1158. start, end = 202, 202
  1159. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1160. assert_almost_equal(fv, fcdyn3[start:end + 1], DECIMAL_4)
  1161. # start nobs, end >nobs 2009q3 - 2015q4
  1162. start, end = 202, 227
  1163. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1164. # start >nobs, end >nobs 2009q4 - 2015q4
  1165. start, end = 203, 227
  1166. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1167. assert_almost_equal(fv, fcdyn4[start:end + 1], DECIMAL_4)
  1168. # defaults
  1169. start, end = None, None
  1170. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1171. assert_almost_equal(fv, fcdyn[5:203], DECIMAL_4)
  1172. def test_arima_predict_css_diffs():
  1173. cpi = load_macrodata_pandas().data['cpi'].values
  1174. # NOTE: Doing no-constant for now to kick the conditional exogenous
  1175. # issue 274 down the road
  1176. # go ahead and git the model to set up necessary variables
  1177. res1 = ARIMA(cpi, (4, 1, 1)).fit(disp=-1, method="css",
  1178. trend="c")
  1179. # but use gretl parameters to predict to avoid precision problems
  1180. params = np.array([0.78349893861244,
  1181. -0.533444105973324,
  1182. 0.321103691668809,
  1183. 0.264012463189186,
  1184. 0.107888256920655,
  1185. 0.920132542916995])
  1186. # we report mean, should we report constant?
  1187. params[0] = params[0] / (1 - params[1:5].sum())
  1188. file_path = os.path.join(current_path, 'results',
  1189. 'results_arima_forecasts_all_css_diff.csv')
  1190. with open(file_path, "rb") as test_data:
  1191. arima_forecasts = np.genfromtxt(test_data, delimiter=",",
  1192. skip_header=1, dtype=float)
  1193. fc = arima_forecasts[:, 0]
  1194. fcdyn = arima_forecasts[:, 1]
  1195. fcdyn2 = arima_forecasts[:, 2]
  1196. fcdyn3 = arima_forecasts[:, 3]
  1197. fcdyn4 = arima_forecasts[:, 4]
  1198. start, end = 1, 3
  1199. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1200. res1.model.predict(params, start, end)
  1201. # start < p, end 0 1959q3 - 1960q1
  1202. start, end = 2, 4
  1203. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1204. res1.model.predict(params, start, end)
  1205. # start < p, end >0 1959q3 - 1971q4
  1206. start, end = 2, 51
  1207. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1208. res1.model.predict(params, start, end)
  1209. # start < p, end nobs 1959q3 - 2009q3
  1210. start, end = 2, 202
  1211. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1212. res1.model.predict(params, start, end)
  1213. # start < p, end >nobs 1959q3 - 2015q4
  1214. start, end = 2, 227
  1215. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1216. res1.model.predict(params, start, end)
  1217. # start 0, end >0 1960q1 - 1971q4
  1218. start, end = 5, 51
  1219. fv = res1.model.predict(params, start, end)
  1220. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1221. # start 0, end nobs 1960q1 - 2009q3
  1222. start, end = 5, 202
  1223. fv = res1.model.predict(params, start, end)
  1224. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1225. # start 0, end >nobs 1960q1 - 2015q4
  1226. # TODO: why detoriating precision?
  1227. fv = res1.model.predict(params, start, end)
  1228. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1229. # start >p, end >0 1965q1 - 1971q4
  1230. start, end = 24, 51
  1231. fv = res1.model.predict(params, start, end)
  1232. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1233. # start >p, end nobs 1965q1 - 2009q3
  1234. start, end = 24, 202
  1235. fv = res1.model.predict(params, start, end)
  1236. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1237. # start >p, end >nobs 1965q1 - 2015q4
  1238. start, end = 24, 227
  1239. fv = res1.model.predict(params, start, end)
  1240. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1241. # start nobs, end nobs 2009q3 - 2009q3
  1242. start, end = 202, 202
  1243. fv = res1.model.predict(params, start, end)
  1244. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1245. # start nobs, end >nobs 2009q3 - 2015q4
  1246. start, end = 202, 227
  1247. fv = res1.model.predict(params, start, end)
  1248. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1249. # start >nobs, end >nobs 2009q4 - 2015q4
  1250. start, end = 203, 227
  1251. fv = res1.model.predict(params, start, end)
  1252. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1253. # defaults
  1254. start, end = None, None
  1255. fv = res1.model.predict(params, start, end)
  1256. assert_almost_equal(fv, fc[5:203], DECIMAL_4)
  1257. # Dynamic
  1258. # start < p, end <p 1959q2 - 1959q4
  1259. start, end = 1, 3
  1260. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1261. res1.predict(start, end, dynamic=True)
  1262. # start < p, end 0 1959q3 - 1960q1
  1263. start, end = 2, 4
  1264. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1265. res1.predict(start, end, dynamic=True)
  1266. # start < p, end >0 1959q3 - 1971q4
  1267. start, end = 2, 51
  1268. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1269. res1.predict(start, end, dynamic=True)
  1270. # start < p, end nobs 1959q3 - 2009q3
  1271. start, end = 2, 202
  1272. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1273. res1.predict(start, end, dynamic=True)
  1274. # start < p, end >nobs 1959q3 - 2015q4
  1275. start, end = 2, 227
  1276. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1277. res1.predict(start, end, dynamic=True)
  1278. # start 0, end >0 1960q1 - 1971q4
  1279. start, end = 5, 51
  1280. fv = res1.model.predict(params, start, end, dynamic=True)
  1281. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1282. # start 0, end nobs 1960q1 - 2009q3
  1283. start, end = 5, 202
  1284. fv = res1.model.predict(params, start, end, dynamic=True)
  1285. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1286. # start 0, end >nobs 1960q1 - 2015q4
  1287. start, end = 5, 227
  1288. fv = res1.model.predict(params, start, end, dynamic=True)
  1289. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1290. # start >p, end >0 1965q1 - 1971q4
  1291. start, end = 24, 51
  1292. fv = res1.model.predict(params, start, end, dynamic=True)
  1293. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1294. # start >p, end nobs 1965q1 - 2009q3
  1295. start, end = 24, 202
  1296. fv = res1.model.predict(params, start, end, dynamic=True)
  1297. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1298. # start >p, end >nobs 1965q1 - 2015q4
  1299. start, end = 24, 227
  1300. fv = res1.model.predict(params, start, end, dynamic=True)
  1301. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1302. # start nobs, end nobs 2009q3 - 2009q3
  1303. start, end = 202, 202
  1304. fv = res1.model.predict(params, start, end, dynamic=True)
  1305. assert_almost_equal(fv, fcdyn3[start:end + 1], DECIMAL_4)
  1306. # start nobs, end >nobs 2009q3 - 2015q4
  1307. start, end = 202, 227
  1308. fv = res1.model.predict(params, start, end, dynamic=True)
  1309. # start >nobs, end >nobs 2009q4 - 2015q4
  1310. start, end = 203, 227
  1311. fv = res1.model.predict(params, start, end, dynamic=True)
  1312. assert_almost_equal(fv, fcdyn4[start:end + 1], DECIMAL_4)
  1313. # defaults
  1314. start, end = None, None
  1315. fv = res1.model.predict(params, start, end, dynamic=True)
  1316. assert_almost_equal(fv, fcdyn[5:203], DECIMAL_4)
  1317. def test_arima_predict_mle_diffs():
  1318. cpi = load_macrodata_pandas().data['cpi'].values
  1319. # NOTE: Doing no-constant for now to kick the conditional exogenous
  1320. # issue 274 down the road
  1321. # go ahead and git the model to set up necessary variables
  1322. res1 = ARIMA(cpi, (4, 1, 1)).fit(disp=-1, trend="c")
  1323. # but use gretl parameters to predict to avoid precision problems
  1324. params = np.array([0.926875951549299,
  1325. -0.555862621524846,
  1326. 0.320865492764400,
  1327. 0.252253019082800,
  1328. 0.113624958031799,
  1329. 0.939144026934634])
  1330. file_path = os.path.join(current_path, 'results',
  1331. 'results_arima_forecasts_all_mle_diff.csv')
  1332. with open(file_path, "rb") as test_data:
  1333. arima_forecasts = np.genfromtxt(test_data, delimiter=",",
  1334. skip_header=1, dtype=float)
  1335. fc = arima_forecasts[:, 0]
  1336. fcdyn = arima_forecasts[:, 1]
  1337. fcdyn2 = arima_forecasts[:, 2]
  1338. fcdyn3 = arima_forecasts[:, 3]
  1339. fcdyn4 = arima_forecasts[:, 4]
  1340. start, end = 1, 3
  1341. res1.model.predict(params, start, end)
  1342. # start < p, end 0 1959q3 - 1960q1
  1343. start, end = 2, 4
  1344. res1.model.predict(params, start, end)
  1345. # start < p, end >0 1959q3 - 1971q4
  1346. start, end = 2, 51
  1347. res1.model.predict(params, start, end)
  1348. # start < p, end nobs 1959q3 - 2009q3
  1349. start, end = 2, 202
  1350. res1.model.predict(params, start, end)
  1351. # start < p, end >nobs 1959q3 - 2015q4
  1352. start, end = 2, 227
  1353. res1.model.predict(params, start, end)
  1354. # start 0, end >0 1960q1 - 1971q4
  1355. start, end = 5, 51
  1356. fv = res1.model.predict(params, start, end)
  1357. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1358. # start 0, end nobs 1960q1 - 2009q3
  1359. start, end = 5, 202
  1360. fv = res1.model.predict(params, start, end)
  1361. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1362. # start 0, end >nobs 1960q1 - 2015q4
  1363. # TODO: why detoriating precision?
  1364. fv = res1.model.predict(params, start, end)
  1365. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1366. # start >p, end >0 1965q1 - 1971q4
  1367. start, end = 24, 51
  1368. fv = res1.model.predict(params, start, end)
  1369. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1370. # start >p, end nobs 1965q1 - 2009q3
  1371. start, end = 24, 202
  1372. fv = res1.model.predict(params, start, end)
  1373. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1374. # start >p, end >nobs 1965q1 - 2015q4
  1375. start, end = 24, 227
  1376. fv = res1.model.predict(params, start, end)
  1377. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1378. # start nobs, end nobs 2009q3 - 2009q3
  1379. start, end = 202, 202
  1380. fv = res1.model.predict(params, start, end)
  1381. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1382. # start nobs, end >nobs 2009q3 - 2015q4
  1383. start, end = 202, 227
  1384. fv = res1.model.predict(params, start, end)
  1385. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1386. # start >nobs, end >nobs 2009q4 - 2015q4
  1387. start, end = 203, 227
  1388. fv = res1.model.predict(params, start, end)
  1389. assert_almost_equal(fv, fc[start:end + 1], DECIMAL_4)
  1390. # defaults
  1391. start, end = None, None
  1392. fv = res1.model.predict(params, start, end)
  1393. assert_almost_equal(fv, fc[1:203], DECIMAL_4)
  1394. # Dynamic
  1395. # NOTE: should raise
  1396. # start < p, end <p 1959q2 - 1959q4
  1397. start, end = 1, 3
  1398. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1399. res1.predict(start, end, dynamic=True)
  1400. # start < p, end 0 1959q3 - 1960q1
  1401. start, end = 2, 4
  1402. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1403. res1.predict(start, end, dynamic=True)
  1404. # start < p, end >0 1959q3 - 1971q4
  1405. start, end = 2, 51
  1406. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1407. res1.predict(start, end, dynamic=True)
  1408. # start < p, end nobs 1959q3 - 2009q3
  1409. start, end = 2, 202
  1410. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1411. res1.predict(start, end, dynamic=True)
  1412. # start < p, end >nobs 1959q3 - 2015q4
  1413. start, end = 2, 227
  1414. with pytest.raises(ValueError, match='Start must be >= k_ar'):
  1415. res1.predict(start, end, dynamic=True)
  1416. # start 0, end >0 1960q1 - 1971q4
  1417. start, end = 5, 51
  1418. fv = res1.model.predict(params, start, end, dynamic=True)
  1419. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1420. # start 0, end nobs 1960q1 - 2009q3
  1421. start, end = 5, 202
  1422. fv = res1.model.predict(params, start, end, dynamic=True)
  1423. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1424. # start 0, end >nobs 1960q1 - 2015q4
  1425. start, end = 5, 227
  1426. fv = res1.model.predict(params, start, end, dynamic=True)
  1427. assert_almost_equal(fv, fcdyn[start:end + 1], DECIMAL_4)
  1428. # start >p, end >0 1965q1 - 1971q4
  1429. start, end = 24, 51
  1430. fv = res1.model.predict(params, start, end, dynamic=True)
  1431. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1432. # start >p, end nobs 1965q1 - 2009q3
  1433. start, end = 24, 202
  1434. fv = res1.model.predict(params, start, end, dynamic=True)
  1435. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1436. # start >p, end >nobs 1965q1 - 2015q4
  1437. start, end = 24, 227
  1438. fv = res1.model.predict(params, start, end, dynamic=True)
  1439. assert_almost_equal(fv, fcdyn2[start:end + 1], DECIMAL_4)
  1440. # start nobs, end nobs 2009q3 - 2009q3
  1441. start, end = 202, 202
  1442. fv = res1.model.predict(params, start, end, dynamic=True)
  1443. assert_almost_equal(fv, fcdyn3[start:end + 1], DECIMAL_4)
  1444. # start nobs, end >nobs 2009q3 - 2015q4
  1445. start, end = 202, 227
  1446. fv = res1.model.predict(params, start, end, dynamic=True)
  1447. # start >nobs, end >nobs 2009q4 - 2015q4
  1448. start, end = 203, 227
  1449. fv = res1.model.predict(params, start, end, dynamic=True)
  1450. assert_almost_equal(fv, fcdyn4[start:end + 1], DECIMAL_4)
  1451. # defaults
  1452. start, end = None, None
  1453. fv = res1.model.predict(params, start, end, dynamic=True)
  1454. assert_almost_equal(fv, fcdyn[5:203], DECIMAL_4)
  1455. def test_arima_wrapper():
  1456. cpi = load_macrodata_pandas().data['cpi']
  1457. cpi.index = pd.Index(cpi_dates)
  1458. res = ARIMA(cpi, (4, 1, 1), freq='Q').fit(disp=-1)
  1459. assert_equal(res.params.index,
  1460. pd.Index(['const', 'ar.L1.D.cpi', 'ar.L2.D.cpi',
  1461. 'ar.L3.D.cpi', 'ar.L4.D.cpi',
  1462. 'ma.L1.D.cpi']))
  1463. assert_equal(res.model.endog_names, 'D.cpi')
  1464. @pytest.mark.smoke
  1465. def test_1dexog():
  1466. # smoke test, this will raise an error if broken
  1467. dta = load_macrodata_pandas().data
  1468. endog = dta['realcons'].values
  1469. exog = dta['m1'].values.squeeze()
  1470. with warnings.catch_warnings():
  1471. warnings.simplefilter("ignore")
  1472. mod = ARMA(endog, (1, 1), exog).fit(disp=-1)
  1473. mod.predict(193, 203, exog[-10:])
  1474. # check for dynamic is true and pandas Series see #2589
  1475. mod.predict(193, 202, exog[-10:], dynamic=True)
  1476. dta.index = pd.Index(cpi_dates)
  1477. mod = ARMA(dta['realcons'], (1, 1), dta['m1']).fit(disp=-1)
  1478. mod.predict(dta.index[-10], dta.index[-1], exog=dta['m1'][-10:],
  1479. dynamic=True)
  1480. mod = ARMA(dta['realcons'], (1, 1), dta['m1']).fit(trend='nc', disp=-1)
  1481. mod.predict(dta.index[-10], dta.index[-1], exog=dta['m1'][-10:],
  1482. dynamic=True)
  1483. def test_arima_predict_bug():
  1484. # predict_start_date was not getting set on start = None
  1485. from statsmodels.datasets import sunspots
  1486. dta = sunspots.load_pandas().data.SUNACTIVITY
  1487. dta.index = pd.date_range(start='1700', end='2009', freq='A')[:309]
  1488. arma_mod20 = ARMA(dta, (2, 0)).fit(disp=-1)
  1489. arma_mod20.predict(None, None)
  1490. # test prediction with time stamp, see #2587
  1491. predict = arma_mod20.predict(dta.index[-20], dta.index[-1])
  1492. assert predict.index.equals(dta.index[-20:])
  1493. predict = arma_mod20.predict(dta.index[-20], dta.index[-1], dynamic=True)
  1494. assert predict.index.equals(dta.index[-20:])
  1495. # partially out of sample
  1496. predict_dates = pd.date_range(start='2000', end='2015', freq='A')
  1497. predict = arma_mod20.predict(predict_dates[0], predict_dates[-1])
  1498. assert predict.index.equals(predict_dates)
  1499. def test_arima_predict_q2():
  1500. # bug with q > 1 for arima predict
  1501. inv = load_macrodata_pandas().data['realinv'].values
  1502. arima_mod = ARIMA(np.log(inv), (1, 1, 2)).fit(start_params=[0, 0, 0, 0],
  1503. disp=-1)
  1504. fc, stderr, conf_int = arima_mod.forecast(5)
  1505. # values copy-pasted from gretl
  1506. assert_almost_equal(fc,
  1507. [7.306320, 7.313825, 7.321749, 7.329827, 7.337962],
  1508. 5)
  1509. def test_arima_predict_pandas_nofreq():
  1510. # this is issue 712
  1511. dates = ["2010-01-04", "2010-01-05", "2010-01-06", "2010-01-07",
  1512. "2010-01-08", "2010-01-11", "2010-01-12", "2010-01-11",
  1513. "2010-01-12", "2010-01-13", "2010-01-17"]
  1514. close = [626.75, 623.99, 608.26, 594.1, 602.02, 601.11, 590.48, 587.09,
  1515. 589.85, 580.0, 587.62]
  1516. data = pd.DataFrame(close, index=DatetimeIndex(dates), columns=["close"])
  1517. # TODO: fix this names bug for non-string names names
  1518. with pytest.warns(ValueWarning, match="it has no associated frequency"):
  1519. arma = ARMA(data, order=(1, 0)).fit(disp=-1)
  1520. # first check that in-sample prediction works
  1521. predict = arma.predict()
  1522. assert predict.index.equals(data.index)
  1523. # check that this raises an exception when date not on index
  1524. assert_raises(KeyError, arma.predict, start="2010-1-9", end=10)
  1525. assert_raises(KeyError, arma.predict, start="2010-1-9", end="2010-1-17")
  1526. # raise because end not on index
  1527. assert_raises(KeyError, arma.predict, start="2010-1-4", end="2010-1-10")
  1528. # raise because end not on index
  1529. assert_raises(KeyError, arma.predict, start=3, end="2010-1-10")
  1530. predict = arma.predict(start="2010-1-7", end=10) # should be of length 10
  1531. assert len(predict) == 8
  1532. assert predict.index.equals(data.index[3:10 + 1])
  1533. with pytest.warns(ValueWarning, match="No supported index is available"):
  1534. predict = arma.predict(start="2010-1-7", end=14)
  1535. assert predict.index.equals(pd.Index(lrange(3, 15)))
  1536. with pytest.warns(ValueWarning, match="No supported index is available"):
  1537. predict = arma.predict(start=3, end=14)
  1538. assert predict.index.equals(pd.Index(lrange(3, 15)))
  1539. # end can be a date if it's in the sample and on the index
  1540. # predict dates is just a slice of the dates index then
  1541. predict = arma.predict(start="2010-1-6", end="2010-1-13")
  1542. assert predict.index.equals(data.index[2:10])
  1543. predict = arma.predict(start=2, end="2010-1-13")
  1544. assert predict.index.equals(data.index[2:10])
  1545. def test_arima_predict_exog():
  1546. # check 625 and 626
  1547. # from statsmodels.tsa.arima_process import arma_generate_sample
  1548. # arparams = np.array([1, -.45, .25])
  1549. # maparams = np.array([1, .15])
  1550. # nobs = 100
  1551. # np.random.seed(123)
  1552. # y = arma_generate_sample(arparams, maparams, nobs, burnin=100)
  1553. # make an exogenous trend
  1554. # X = np.array(lrange(nobs)) / 20.0
  1555. # add a constant
  1556. # y += 2.5
  1557. from pandas import read_csv
  1558. file_name = "/results/results_arima_exog_forecasts_mle.csv"
  1559. arima_forecasts = read_csv(current_path + file_name)
  1560. y = arima_forecasts["y"].dropna()
  1561. X = np.arange(len(y) + 25) / 20.
  1562. predict_expected = arima_forecasts["predict"]
  1563. arma_res = ARMA(y.values, order=(2, 1), exog=X[:100]).fit(trend="c",
  1564. disp=-1)
  1565. # params from gretl
  1566. params = np.array([2.786912485145725, -0.122650190196475,
  1567. 0.533223846028938, -0.319344321763337,
  1568. 0.132883233000064])
  1569. assert_almost_equal(arma_res.params, params, 5)
  1570. # no exog for in-sample
  1571. predict = arma_res.predict()
  1572. assert_almost_equal(predict, predict_expected.values[:100], 5)
  1573. # check 626
  1574. assert len(arma_res.model.exog_names) == 5
  1575. # exog for out-of-sample and in-sample dynamic
  1576. predict = arma_res.model.predict(params, end=124, exog=X[100:])
  1577. assert_almost_equal(predict, predict_expected.values, 6)
  1578. # conditional sum of squares
  1579. # arima_forecasts = read_csv(current_path + "/results/"
  1580. # "results_arima_exog_forecasts_css.csv")
  1581. # predict_expected = arima_forecasts["predict"].dropna()
  1582. # arma_res = ARMA(y.values, order=(2,1), exog=X[:100]).fit(trend="c",
  1583. # method="css",
  1584. # disp=-1)
  1585. # params = np.array([2.152350033809826, -0.103602399018814,
  1586. # 0.566716580421188, -0.326208009247944,
  1587. # 0.102142932143421])
  1588. # predict = arma_res.model.predict(params)
  1589. # in-sample
  1590. # assert_almost_equal(predict, predict_expected.values[:98], 6)
  1591. # predict = arma_res.model.predict(params, end=124, exog=X[100:])
  1592. # exog for out-of-sample and in-sample dynamic
  1593. # assert_almost_equal(predict, predict_expected.values, 3)
  1594. @pytest.mark.smoke
  1595. def test_arima_no_diff():
  1596. # GH#736
  1597. # smoke test, predict will break if we have ARIMAResults but
  1598. # ARMA model, need ARIMA(p, 0, q) to return an ARMA in init.
  1599. ar = [1, -.75, .15, .35]
  1600. ma = [1, .25, .9]
  1601. np.random.seed(12345)
  1602. y = arma_generate_sample(ar, ma, 100)
  1603. mod = ARIMA(y, (3, 0, 2))
  1604. assert type(mod) is ARMA
  1605. res = mod.fit(disp=-1)
  1606. # smoke test just to be sure
  1607. res.predict()
  1608. @pytest.mark.smoke
  1609. def test_arima_predict_noma():
  1610. # GH#657
  1611. ar = [1, .75]
  1612. ma = [1]
  1613. np.random.seed(12345)
  1614. data = arma_generate_sample(ar, ma, 100)
  1615. arma = ARMA(data, order=(0, 1))
  1616. arma_res = arma.fit(disp=-1)
  1617. arma_res.forecast(1)
  1618. def test_arimax():
  1619. dta = load_macrodata_pandas().data
  1620. dta.index = cpi_dates
  1621. dta = dta[["realdpi", "m1", "realgdp"]]
  1622. y = dta.pop("realdpi")
  1623. # 1 exog
  1624. # X = dta.iloc[1:]["m1"]
  1625. # res = ARIMA(y, (2, 1, 1), X).fit(disp=-1)
  1626. # params = [23.902305009084373, 0.024650911502790, -0.162140641341602,
  1627. # 0.165262136028113, -0.066667022903974]
  1628. # assert_almost_equal(res.params.values, params, 6)
  1629. # 2 exog
  1630. X = dta
  1631. res = ARIMA(y, (2, 1, 1), X).fit(disp=False, solver="nm", maxiter=1000,
  1632. ftol=1e-12, xtol=1e-12)
  1633. # from gretl
  1634. # params = [13.113976653926638, -0.003792125069387, 0.004123504809217,
  1635. # -0.199213760940898, 0.151563643588008, -0.033088661096699]
  1636. # from stata using double
  1637. stata_llf = -1076.108614859121
  1638. params = [13.1259220104, -0.00376814509403812, 0.00411970083135622,
  1639. -0.19921477896158524, 0.15154396192855729, -0.03308400760360837]
  1640. # we can get close
  1641. assert_almost_equal(res.params.values, params, 4)
  1642. # This shows that it's an optimizer problem and not a problem in the code
  1643. assert_almost_equal(res.model.loglike(np.array(params)), stata_llf, 6)
  1644. X = dta.diff()
  1645. X.iloc[0] = 0
  1646. res = ARIMA(y, (2, 1, 1), X).fit(disp=False)
  1647. # gretl will not estimate this - looks like maybe a bug on their part,
  1648. # but we can just fine, we're close to Stata's answer
  1649. # from Stata
  1650. params = [19.5656863783347, 0.32653841355833396198,
  1651. 0.36286527042965188716, -1.01133792126884,
  1652. -0.15722368379307766206, 0.69359822544092153418]
  1653. assert_almost_equal(res.params.values, params, 3)
  1654. def test_bad_start_params():
  1655. endog = np.array([
  1656. 820.69093, 781.0103028, 785.8786988, 767.64282267,
  1657. 778.9837648, 824.6595702, 813.01877867, 751.65598567,
  1658. 753.431091, 746.920813, 795.6201904, 772.65732833,
  1659. 793.4486454, 868.8457766, 823.07226547, 783.09067747,
  1660. 791.50723847, 770.93086347, 835.34157333, 810.64147947,
  1661. 738.36071367, 776.49038513, 822.93272333, 815.26461227,
  1662. 773.70552987, 777.3726522, 811.83444853, 840.95489133,
  1663. 777.51031933, 745.90077307, 806.95113093, 805.77521973,
  1664. 756.70927733, 749.89091773, 1694.2266924, 2398.4802244,
  1665. 1434.6728516, 909.73940427, 929.01291907, 769.07561453,
  1666. 801.1112548, 796.16163313, 817.2496376, 857.73046447,
  1667. 838.849345, 761.92338873, 731.7842242, 770.4641844])
  1668. mod = ARMA(endog, (15, 0))
  1669. with pytest.raises(ValueError, match="pass your own start_params"):
  1670. with pytest.warns(RuntimeWarning, match="encountered in sqrt"):
  1671. # numpy warns "invalid value encountered in sqrt"
  1672. mod.fit()
  1673. inv = load_macrodata_pandas().data['realinv'].values
  1674. arima_mod = ARIMA(np.log(inv), (1, 1, 2))
  1675. with pytest.raises(ValueError):
  1676. arima_mod.fit()
  1677. def test_arima_small_data_bug():
  1678. # Issue 1038, too few observations with given order
  1679. import statsmodels.api as sm
  1680. vals = [96.2, 98.3, 99.1, 95.5, 94.0, 87.1, 87.9, 86.7402777504474]
  1681. dr = pd.date_range("1990", periods=len(vals), freq='Q')
  1682. ts = pd.Series(vals, index=dr)
  1683. df = pd.DataFrame(ts)
  1684. mod = sm.tsa.ARIMA(df, (2, 0, 2))
  1685. assert_raises(ValueError, mod.fit)
  1686. @pytest.mark.smoke
  1687. def test_arima_dataframe_integer_name():
  1688. # Smoke Test for Issue GH#1038
  1689. import statsmodels.api as sm
  1690. vals = [96.2, 98.3, 99.1, 95.5, 94.0, 87.1, 87.9, 86.7402777504474,
  1691. 94.0, 96.5, 93.3, 97.5, 96.3, 92.]
  1692. dr = pd.date_range("1990", periods=len(vals), freq='Q')
  1693. ts = pd.Series(vals, index=dr)
  1694. df = pd.DataFrame(ts)
  1695. sm.tsa.ARIMA(df, (2, 0, 2))
  1696. def test_arima_exog_predict_1d():
  1697. # test 1067
  1698. np.random.seed(12345)
  1699. y = np.random.random(100)
  1700. x = np.random.random(100)
  1701. mod = ARMA(y, (2, 1), x).fit(disp=-1)
  1702. newx = np.random.random(10)
  1703. mod.forecast(steps=10, alpha=0.05, exog=newx)
  1704. with pytest.raises(ValueError):
  1705. mod.forecast(steps=10, alpha=0.05, exog=newx[:5])
  1706. with pytest.raises(ValueError):
  1707. mod.forecast(steps=10, alpha=0.05)
  1708. too_many = pd.DataFrame(np.zeros((10, 2)),
  1709. columns=['x1', 'x2'])
  1710. with pytest.raises(ValueError):
  1711. mod.forecast(steps=10, alpha=0.05, exog=too_many)
  1712. def test_arima_1123():
  1713. # test ARMAX predict when trend is none
  1714. np.random.seed(12345)
  1715. arparams = np.array([.75, -.25])
  1716. maparams = np.array([.65, .35])
  1717. nobs = 20
  1718. np.random.seed(12345)
  1719. y = arma_generate_sample(arparams, maparams, nobs)
  1720. X = np.random.randn(nobs)
  1721. y += 5 * X
  1722. mod = ARMA(y[:-1], order=(1, 0), exog=X[:-1])
  1723. res = mod.fit(trend='nc', disp=False)
  1724. fc = res.forecast(exog=X[-1:])
  1725. # results from gretl
  1726. assert_almost_equal(fc[0], 2.200393, 6)
  1727. assert_almost_equal(fc[1], 1.030743, 6)
  1728. assert_almost_equal(fc[2][0, 0], 0.180175, 6)
  1729. assert_almost_equal(fc[2][0, 1], 4.220611, 6)
  1730. mod = ARMA(y[:-1], order=(1, 1), exog=X[:-1])
  1731. res = mod.fit(trend='nc', disp=False)
  1732. fc = res.forecast(exog=X[-1:])
  1733. assert_almost_equal(fc[0], 2.765688, 6)
  1734. assert_almost_equal(fc[1], 0.835048, 6)
  1735. assert_almost_equal(fc[2][0, 0], 1.129023, 6)
  1736. assert_almost_equal(fc[2][0, 1], 4.402353, 6)
  1737. # make sure this works to. code looked fishy.
  1738. mod = ARMA(y[:-1], order=(1, 0), exog=X[:-1])
  1739. res = mod.fit(trend='c', disp=False)
  1740. fc = res.forecast(exog=X[-1:])
  1741. assert_almost_equal(fc[0], 2.481219, 6)
  1742. assert_almost_equal(fc[1], 0.968759, 6)
  1743. assert_almost_equal(fc[2][0], [0.582485, 4.379952], 6)
  1744. def test_small_data():
  1745. # 1146
  1746. y = [-1214.360173, -1848.209905, -2100.918158, -3647.483678, -4711.186773]
  1747. # refuse to estimate these
  1748. assert_raises(ValueError, ARIMA, y, (2, 0, 3))
  1749. assert_raises(ValueError, ARIMA, y, (1, 1, 3))
  1750. mod = ARIMA(y, (1, 0, 3))
  1751. assert_raises(ValueError, mod.fit, trend="c")
  1752. # try to estimate these...leave it up to the user to check for garbage
  1753. # and be clear, these are garbage parameters.
  1754. # X-12 arima will estimate, gretl refuses to estimate likely a problem
  1755. # in start params regression.
  1756. mod.fit(trend="nc", disp=0, start_params=[.1, .1, .1, .1])
  1757. mod = ARIMA(y, (1, 0, 2))
  1758. import warnings
  1759. with warnings.catch_warnings():
  1760. warnings.simplefilter("ignore")
  1761. mod.fit(disp=0, start_params=[np.mean(y), .1, .1, .1])
  1762. class TestARMA00(object):
  1763. @classmethod
  1764. def setup_class(cls):
  1765. from statsmodels.datasets.sunspots import load_pandas
  1766. sunspots = load_pandas().data['SUNACTIVITY'].values
  1767. cls.y = y = sunspots
  1768. cls.arma_00_model = ARMA(y, order=(0, 0))
  1769. cls.arma_00_res = cls.arma_00_model.fit(disp=-1)
  1770. def test_parameters(self):
  1771. params = self.arma_00_res.params
  1772. assert_almost_equal(self.y.mean(), params)
  1773. def test_predictions(self):
  1774. predictions = self.arma_00_res.predict()
  1775. assert_almost_equal(self.y.mean() * np.ones_like(predictions),
  1776. predictions)
  1777. def test_arroots(self):
  1778. # regression test; older implementation of arroots returned None
  1779. # instead of en empty array
  1780. roots = self.arma_00_res.arroots
  1781. assert_equal(roots.size, 0)
  1782. def test_maroots(self):
  1783. # regression test; older implementation of arroots returned None
  1784. # instead of en empty array
  1785. roots = self.arma_00_res.maroots
  1786. assert_equal(roots.size, 0)
  1787. @pytest.mark.skip("This test is invalid since the ICs differ due to "
  1788. "df_model differences between OLS and ARIMA")
  1789. def test_information_criteria(self):
  1790. res = self.arma_00_res
  1791. y = self.y
  1792. ols_res = OLS(y, np.ones_like(y)).fit(disp=-1)
  1793. ols_ic = np.array([ols_res.aic, ols_res.bic])
  1794. arma_ic = np.array([res.aic, res.bic])
  1795. assert_almost_equal(ols_ic, arma_ic, DECIMAL_4)
  1796. def test_arma_00_nc(self):
  1797. arma_00 = ARMA(self.y, order=(0, 0))
  1798. assert_raises(ValueError, arma_00.fit, trend='nc', disp=-1)
  1799. def test_css(self):
  1800. arma = ARMA(self.y, order=(0, 0))
  1801. fit = arma.fit(method='css', disp=-1)
  1802. predictions = fit.predict()
  1803. assert_almost_equal(self.y.mean() * np.ones_like(predictions),
  1804. predictions)
  1805. def test_arima(self):
  1806. yi = np.cumsum(self.y)
  1807. arima = ARIMA(yi, order=(0, 1, 0))
  1808. fit = arima.fit(disp=-1)
  1809. assert_almost_equal(np.diff(yi).mean(), fit.params, DECIMAL_4)
  1810. def test_arma_ols(self):
  1811. y = self.y
  1812. y_lead = y[1:]
  1813. y_lag = y[:-1]
  1814. T = y_lag.shape[0]
  1815. X = np.hstack((np.ones((T, 1)), y_lag[:, None]))
  1816. ols_res = OLS(y_lead, X).fit()
  1817. arma_res = ARMA(y_lead, order=(0, 0), exog=y_lag).fit(trend='c',
  1818. disp=-1)
  1819. assert_almost_equal(ols_res.params, arma_res.params)
  1820. def test_arma_exog_no_constant(self):
  1821. y = self.y
  1822. y_lead = y[1:]
  1823. y_lag = y[:-1]
  1824. X = y_lag[:, None]
  1825. ols_res = OLS(y_lead, X).fit()
  1826. arma_res = ARMA(y_lead, order=(0, 0), exog=y_lag).fit(trend='nc',
  1827. disp=-1)
  1828. assert_almost_equal(ols_res.params, arma_res.params)
  1829. def test_arima_dates_startatend():
  1830. # bug
  1831. np.random.seed(18)
  1832. x = pd.Series(np.random.random(36),
  1833. index=pd.date_range(start='1/1/1990',
  1834. periods=36, freq='M'))
  1835. res = ARIMA(x, (1, 0, 0)).fit(disp=0)
  1836. pred = res.predict(start=len(x), end=len(x))
  1837. assert pred.index[0] == x.index.shift(1)[-1]
  1838. fc = res.forecast()[0]
  1839. assert_almost_equal(pred.values[0], fc)
  1840. def test_arma_missing():
  1841. from statsmodels.tools.sm_exceptions import MissingDataError
  1842. # bug 1343
  1843. y = np.random.random(40)
  1844. y[-1] = np.nan
  1845. assert_raises(MissingDataError, ARMA, y, (1, 0), missing='raise')
  1846. @pytest.mark.matplotlib
  1847. def test_plot_predict(close_figures):
  1848. from statsmodels.datasets.sunspots import load_pandas
  1849. dta = load_pandas().data[['SUNACTIVITY']]
  1850. dta.index = date_range(start='1700', end='2009', freq='A')[:309]
  1851. res = ARMA(dta, (3, 0)).fit(disp=-1)
  1852. res.plot_predict('1990', '2012', dynamic=True, plot_insample=False)
  1853. res.plot_predict('1990', '2012', dynamic=True, plot_insample=True)
  1854. res = ARIMA(dta, (3, 1, 0)).fit(disp=-1)
  1855. res.plot_predict('1990', '2012', dynamic=True, plot_insample=False)
  1856. res.plot_predict('1990', '2012', dynamic=True, plot_insample=True)
  1857. def test_arima_diff2():
  1858. dta = load_macrodata_pandas().data['cpi']
  1859. dta.index = cpi_dates
  1860. mod = ARIMA(dta, (3, 2, 1)).fit(disp=-1)
  1861. fc, fcerr, conf_int = mod.forecast(10)
  1862. # forecasts from gretl
  1863. conf_int_res = [(216.139, 219.231),
  1864. (216.472, 221.520),
  1865. (217.064, 223.649),
  1866. (217.586, 225.727),
  1867. (218.119, 227.770),
  1868. (218.703, 229.784),
  1869. (219.306, 231.777),
  1870. (219.924, 233.759),
  1871. (220.559, 235.735),
  1872. (221.206, 237.709)]
  1873. fc_res = [217.685, 218.996, 220.356, 221.656, 222.945, 224.243, 225.541,
  1874. 226.841, 228.147, 229.457]
  1875. fcerr_res = [0.7888, 1.2878, 1.6798, 2.0768, 2.4620, 2.8269, 3.1816,
  1876. 3.52950, 3.8715, 4.2099]
  1877. assert_almost_equal(fc, fc_res, 3)
  1878. assert_almost_equal(fcerr, fcerr_res, 3)
  1879. assert_almost_equal(conf_int, conf_int_res, 3)
  1880. predicted = mod.predict('2008Q1', '2012Q1', typ='levels')
  1881. predicted_res = [214.464, 215.478, 221.277, 217.453, 212.419, 213.530,
  1882. 215.087, 217.685, 218.996, 220.356, 221.656,
  1883. 222.945, 224.243, 225.541, 226.841, 228.147,
  1884. 229.457]
  1885. assert_almost_equal(predicted, predicted_res, 3)
  1886. def test_arima111_predict_exog_2127():
  1887. # regression test for issue #2127
  1888. ef = [0.03005, 0.03917, 0.02828, 0.03644, 0.03379, 0.02744,
  1889. 0.03343, 0.02621, 0.03050, 0.02455, 0.03261, 0.03507,
  1890. 0.02734, 0.05373, 0.02677, 0.03443, 0.03331, 0.02741,
  1891. 0.03709, 0.02113, 0.03343, 0.02011, 0.03675, 0.03077,
  1892. 0.02201, 0.04844, 0.05518, 0.03765, 0.05433, 0.03049,
  1893. 0.04829, 0.02936, 0.04421, 0.02457, 0.04007, 0.03009,
  1894. 0.04504, 0.05041, 0.03651, 0.02719, 0.04383, 0.02887,
  1895. 0.03440, 0.03348, 0.02364, 0.03496, 0.02549, 0.03284,
  1896. 0.03523, 0.02579, 0.03080, 0.01784, 0.03237, 0.02078,
  1897. 0.03508, 0.03062, 0.02006, 0.02341, 0.02223, 0.03145,
  1898. 0.03081, 0.02520, 0.02683, 0.01720, 0.02225, 0.01579,
  1899. 0.02237, 0.02295, 0.01830, 0.02356, 0.02051, 0.02932,
  1900. 0.03025, 0.02390, 0.02635, 0.01863, 0.02994, 0.01762,
  1901. 0.02837, 0.02421, 0.01951, 0.02149, 0.02079, 0.02528,
  1902. 0.02575, 0.01634, 0.02563, 0.01719, 0.02915, 0.01724,
  1903. 0.02804, 0.02750, 0.02099, 0.02522, 0.02422, 0.03254,
  1904. 0.02095, 0.03241, 0.01867, 0.03998, 0.02212, 0.03034,
  1905. 0.03419, 0.01866, 0.02623, 0.02052]
  1906. ue = [4.9, 5.0, 5.0, 5.0, 4.9, 4.7, 4.8, 4.7, 4.7,
  1907. 4.6, 4.6, 4.7, 4.7, 4.5, 4.4, 4.5, 4.4, 4.6,
  1908. 4.5, 4.4, 4.5, 4.4, 4.6, 4.7, 4.6, 4.7, 4.7,
  1909. 4.7, 5.0, 5.0, 4.9, 5.1, 5.0, 5.4, 5.6, 5.8,
  1910. 6.1, 6.1, 6.5, 6.8, 7.3, 7.8, 8.3, 8.7, 9.0,
  1911. 9.4, 9.5, 9.5, 9.6, 9.8, 10., 9.9, 9.9, 9.7,
  1912. 9.8, 9.9, 9.9, 9.6, 9.4, 9.5, 9.5, 9.5, 9.5,
  1913. 9.8, 9.4, 9.1, 9.0, 9.0, 9.1, 9.0, 9.1, 9.0,
  1914. 9.0, 9.0, 8.8, 8.6, 8.5, 8.2, 8.3, 8.2, 8.2,
  1915. 8.2, 8.2, 8.2, 8.1, 7.8, 7.8, 7.8, 7.9, 7.9,
  1916. 7.7, 7.5, 7.5, 7.5, 7.5, 7.3, 7.2, 7.2, 7.2,
  1917. 7.0, 6.7, 6.6, 6.7, 6.7, 6.3, 6.3]
  1918. ue = np.array(ue) / 100
  1919. model = ARIMA(ef, (1, 1, 1), exog=ue)
  1920. res = model.fit(transparams=False, pgtol=1e-8, iprint=0, disp=0)
  1921. assert_equal(res.mle_retvals['warnflag'], 0)
  1922. predicts = res.predict(start=len(ef), end=len(ef) + 10,
  1923. exog=ue[-11:], typ='levels')
  1924. # regression test, not verified numbers
  1925. predicts_res = np.array([
  1926. 0.02591095, 0.02321325, 0.02436579, 0.02368759, 0.02389753,
  1927. 0.02372, 0.0237481, 0.0236738, 0.023644, 0.0236283,
  1928. 0.02362267])
  1929. assert_allclose(predicts, predicts_res, atol=5e-6)
  1930. # Smoke check of forecast with exog in ARIMA
  1931. res.forecast(steps=10, exog=np.empty(10))
  1932. with pytest.raises(ValueError):
  1933. res.forecast(steps=10)
  1934. with pytest.raises(ValueError):
  1935. res.forecast(steps=10, exog=np.empty((10, 2)))
  1936. with pytest.raises(ValueError):
  1937. res.forecast(steps=10, exog=np.empty(100))
  1938. def test_arima_exog_predict():
  1939. # test forecasting and dynamic prediction with exog against Stata
  1940. dta = load_macrodata_pandas().data
  1941. cpi_dates = period_range(start='1959Q1', end='2009Q3', freq='Q')
  1942. dta.index = cpi_dates
  1943. data = dta
  1944. data['loginv'] = np.log(data['realinv'])
  1945. data['loggdp'] = np.log(data['realgdp'])
  1946. data['logcons'] = np.log(data['realcons'])
  1947. forecast_period = period_range(start='2008Q2', end='2009Q4', freq='Q')
  1948. end = forecast_period[0]
  1949. data_sample = data.loc[dta.index < end]
  1950. exog_full = data[['loggdp', 'logcons']]
  1951. # pandas
  1952. mod = ARIMA(data_sample['loginv'], (1, 0, 1),
  1953. exog=data_sample[['loggdp', 'logcons']])
  1954. with warnings.catch_warnings():
  1955. warnings.simplefilter("ignore")
  1956. res = mod.fit(disp=0, solver='bfgs', maxiter=5000)
  1957. predicted_arma_fp = res.predict(start=197, end=202,
  1958. exog=exog_full.values[197:]).values
  1959. predicted_arma_dp = res.predict(start=193, end=202,
  1960. exog=exog_full[197:], dynamic=True)
  1961. # numpy
  1962. mod2 = ARIMA(np.asarray(data_sample['loginv']), (1, 0, 1),
  1963. exog=np.asarray(data_sample[['loggdp', 'logcons']]))
  1964. res2 = mod2.fit(start_params=res.params, disp=0,
  1965. solver='bfgs', maxiter=5000)
  1966. exog_full = data[['loggdp', 'logcons']]
  1967. predicted_arma_f = res2.predict(start=197, end=202,
  1968. exog=exog_full.values[197:])
  1969. predicted_arma_d = res2.predict(start=193, end=202,
  1970. exog=exog_full[197:], dynamic=True)
  1971. endog_scale = 100
  1972. ex_scale = 1000
  1973. # ARIMA(1, 1, 1)
  1974. ex = ex_scale * np.asarray(data_sample[['loggdp', 'logcons']].diff())
  1975. # The first observation is not (supposed to be) used,
  1976. # but I get a Lapack problem
  1977. # Intel MKL ERROR: Parameter 5 was incorrect on entry to DLASCL.
  1978. ex[0] = 0
  1979. mod111 = ARIMA(100 * np.asarray(data_sample['loginv']), (1, 1, 1),
  1980. # Stata differences also the exog
  1981. exog=ex)
  1982. sv = [-0.94771574, 0.5869353, -0.32651653, 0.78066562, -0.68747501]
  1983. res111 = mod111.fit(start_params=sv, disp=0, solver='bfgs', maxiter=5000)
  1984. exog_full_d = ex_scale * data[['loggdp', 'logcons']].diff()
  1985. res111.predict(start=197, end=202, exog=exog_full_d.values[197:])
  1986. predicted_arima_f = res111.predict(start=196, end=202,
  1987. exog=exog_full_d.values[197:],
  1988. typ='levels')
  1989. predicted_arima_d = res111.predict(start=193, end=202,
  1990. exog=exog_full_d.values[197:],
  1991. typ='levels', dynamic=True)
  1992. res_f101 = np.array(
  1993. [7.73975859954, 7.71660108543, 7.69808978329, 7.70872117504,
  1994. 7.65183927580, 7.69784279784, 7.70290907856, 7.69237782644,
  1995. 7.65017785174, 7.66061689028, 7.65980022857, 7.61505314129,
  1996. 7.51697158428, 7.51657606630, 7.5271053284])
  1997. res_f111 = np.array(
  1998. [7.74460013693, 7.71958207517, 7.69629561172, 7.71208186737,
  1999. 7.65758850178, 7.69223472572, 7.70411775588, 7.68896109499,
  2000. 7.64016249001, 7.64871881901, 7.62550283402, 7.55814609462,
  2001. 7.44431310053, 7.42963968062, 7.43554675427])
  2002. res_d111 = np.array(
  2003. [7.74460013693, 7.71958207517, 7.69629561172, 7.71208186737,
  2004. 7.65758850178, 7.69223472572, 7.71870821151, 7.72994302150,
  2005. 7.71439447355, 7.72544001101, 7.70521902623, 7.64020040524,
  2006. 7.52819271910, 7.51494426940, 7.52196378005])
  2007. res_d101 = np.array(
  2008. [7.73975859954, 7.71660108543, 7.69808978329, 7.70872117504,
  2009. 7.65183927580, 7.69784279784, 7.72522142662, 7.73962377858,
  2010. 7.73245950636, 7.74935432862, 7.74449584691, 7.69589103679,
  2011. 7.59412746880, 7.59021764836, 7.59739267775])
  2012. # Relax tolerance on OSX due to frequent small failures
  2013. # GH 4657
  2014. tol = 1e-3 if PLATFORM_OSX or PLATFORM_WIN else 1e-4
  2015. assert_allclose(predicted_arma_dp,
  2016. res_d101[-len(predicted_arma_d):], rtol=tol, atol=tol)
  2017. assert_allclose(predicted_arma_fp,
  2018. res_f101[-len(predicted_arma_f):], rtol=tol, atol=tol)
  2019. assert_allclose(predicted_arma_d,
  2020. res_d101[-len(predicted_arma_d):], rtol=tol, atol=tol)
  2021. assert_allclose(predicted_arma_f,
  2022. res_f101[-len(predicted_arma_f):], rtol=tol, atol=tol)
  2023. assert_allclose(predicted_arima_d / endog_scale,
  2024. res_d111[-len(predicted_arima_d):], rtol=tol, atol=tol)
  2025. assert_allclose(predicted_arima_f / endog_scale,
  2026. res_f111[-len(predicted_arima_f):], rtol=tol, atol=tol)
  2027. # test for forecast with 0 ar fix in #2457 numbers again from Stata
  2028. res_f002 = np.array(
  2029. [7.70178181209, 7.67445481224, 7.67153737650, 7.67729153190,
  2030. 7.61173201163, 7.67913499878, 7.67276092120, 7.66275451925,
  2031. 7.65199799315, 7.65149983741, 7.65554131408, 7.62213286298,
  2032. 7.53795983357, 7.53626130154, 7.54539963934])
  2033. res_d002 = np.array(
  2034. [7.70178181209, 7.67445481224, 7.67153737650, 7.67729153190,
  2035. 7.61173201163, 7.67913499878, 7.67306697759, 7.65287924998,
  2036. 7.64904451605, 7.66580449603, 7.66252081172, 7.62213286298,
  2037. 7.53795983357, 7.53626130154, 7.54539963934])
  2038. mod_002 = ARIMA(np.asarray(data_sample['loginv']), (0, 0, 2),
  2039. exog=np.asarray(data_sample[['loggdp', 'logcons']]))
  2040. # does not converge with default starting values
  2041. start_params = np.concatenate((res.params[[0, 1, 2, 4]], [0]))
  2042. res_002 = mod_002.fit(start_params=start_params, disp=0,
  2043. solver='bfgs', maxiter=5000)
  2044. # forecast
  2045. fpredict_002 = res_002.predict(start=197, end=202,
  2046. exog=exog_full.values[197:])
  2047. forecast_002 = res_002.forecast(steps=len(exog_full.values[197:]),
  2048. exog=exog_full.values[197:])
  2049. # TODO: Checking other results
  2050. forecast_002 = forecast_002[0]
  2051. assert_allclose(fpredict_002, res_f002[-len(fpredict_002):],
  2052. rtol=tol, atol=tol / 100.0)
  2053. assert_allclose(forecast_002, res_f002[-len(forecast_002):],
  2054. rtol=tol, atol=tol / 100.0)
  2055. # dynamic predict
  2056. dpredict_002 = res_002.predict(start=193, end=202,
  2057. exog=exog_full.values[197:], dynamic=True)
  2058. assert_allclose(dpredict_002, res_d002[-len(dpredict_002):],
  2059. rtol=tol, atol=tol / 100.0)
  2060. # in-sample dynamic predict should not need exog, #2982
  2061. predict_3a = res_002.predict(start=100, end=120, dynamic=True)
  2062. predict_3b = res_002.predict(start=100, end=120,
  2063. exog=exog_full.values[100:120], dynamic=True)
  2064. assert_allclose(predict_3a, predict_3b, rtol=1e-10)
  2065. h = len(exog_full.values[197:])
  2066. with pytest.raises(ValueError):
  2067. res_002.forecast(steps=h)
  2068. with pytest.raises(ValueError):
  2069. res_002.forecast(steps=h, exog=np.empty((h, 20)))
  2070. with pytest.raises(ValueError):
  2071. res_002.forecast(steps=h, exog=np.empty(20))
  2072. def test_arima_fit_multiple_calls():
  2073. y = [-1214.360173, -1848.209905, -2100.918158, -3647.483678, -4711.186773]
  2074. mod = ARIMA(y, (1, 0, 2))
  2075. # Make multiple calls to fit
  2076. with pytest.warns(HessianInversionWarning, match="no bse or cov"):
  2077. mod.fit(disp=0, start_params=[np.mean(y), .1, .1, .1])
  2078. assert_equal(mod.exog_names, ['const', 'ar.L1.y', 'ma.L1.y', 'ma.L2.y'])
  2079. with pytest.warns(HessianInversionWarning, match="no bse or cov"):
  2080. res = mod.fit(disp=0, start_params=[np.mean(y), .1, .1, .1])
  2081. assert_equal(mod.exog_names, ['const', 'ar.L1.y', 'ma.L1.y', 'ma.L2.y'])
  2082. # ensure summary() works # TODO: separate and pytest.mark.smoke
  2083. res.summary()
  2084. res.summary2()
  2085. # test multiple calls when there is only a constant term
  2086. mod = ARIMA(y, (0, 0, 0))
  2087. # Make multiple calls to fit
  2088. mod.fit(disp=0, start_params=[np.mean(y)])
  2089. assert_equal(mod.exog_names, ['const'])
  2090. with pytest.raises(RuntimeError):
  2091. mod.fit(disp=0, start_params=[np.mean(y)])
  2092. def test_long_ar_start_params():
  2093. np.random.seed(12345)
  2094. arparams = np.array([1, -.75, .25])
  2095. maparams = np.array([1, .65, .35])
  2096. nobs = 30
  2097. np.random.seed(12345)
  2098. y = arma_generate_sample(arparams, maparams, nobs)
  2099. model = ARMA(y, order=(2, 2))
  2100. model.fit(method='css', start_ar_lags=10, disp=0)
  2101. with pytest.raises(RuntimeError):
  2102. model.fit(method='css-mle', start_ar_lags=10, disp=0)
  2103. with pytest.raises(RuntimeError):
  2104. model.fit(method='mle', start_ar_lags=10, disp=0)
  2105. with pytest.raises(RuntimeError):
  2106. model.fit(ValueError, start_ar_lags=nobs + 5, disp=0)
  2107. def test_start_params_too_large():
  2108. np.random.seed(12345)
  2109. arparams = np.array([1, -.75, .25])
  2110. maparams = np.array([1, .65, .35])
  2111. nobs = 30
  2112. np.random.seed(12345)
  2113. y = arma_generate_sample(arparams, maparams, nobs)
  2114. model = ARMA(y, order=(2, 2))
  2115. with pytest.raises(ValueError):
  2116. model.fit(start_ar_lags=nobs + 5, disp=0)
  2117. def test_arma_pickle():
  2118. np.random.seed(9876565)
  2119. x = fa.ArmaFft([1, -0.5], [1., 0.4], 40).generate_sample(nsample=200,
  2120. burnin=1000)
  2121. mod = ARMA(x, (1, 1))
  2122. pkl_mod = pickle.loads(pickle.dumps(mod))
  2123. res = mod.fit(trend="c", disp=-1, solver='newton')
  2124. pkl_res = pkl_mod.fit(trend="c", disp=-1, solver='newton')
  2125. assert_allclose(res.params, pkl_res.params)
  2126. assert_allclose(res.llf, pkl_res.llf)
  2127. assert_almost_equal(res.resid, pkl_res.resid)
  2128. assert_almost_equal(res.fittedvalues, pkl_res.fittedvalues)
  2129. assert_almost_equal(res.pvalues, pkl_res.pvalues)
  2130. def test_arima_pickle():
  2131. endog = y_arma[:, 6]
  2132. mod = ARIMA(endog, (1, 1, 1))
  2133. pkl_mod = pickle.loads(pickle.dumps(mod))
  2134. res = mod.fit(trend="c", disp=-1, solver='newton')
  2135. pkl_res = pkl_mod.fit(trend="c", disp=-1, solver='newton')
  2136. assert_allclose(res.params, pkl_res.params)
  2137. assert_allclose(res.llf, pkl_res.llf)
  2138. assert_almost_equal(res.resid, pkl_res.resid)
  2139. assert_almost_equal(res.fittedvalues, pkl_res.fittedvalues)
  2140. assert_almost_equal(res.pvalues, pkl_res.pvalues)
  2141. def test_arima_not_implemented():
  2142. formula = ' WUE ~ 1 + SFO3 '
  2143. data = [-1214.360173, -1848.209905, -2100.918158]
  2144. assert_raises(NotImplementedError, ARIMA.from_formula, formula, data)
  2145. def test_endog_int():
  2146. # int endog should produce same result as float, #3504
  2147. np.random.seed(123987)
  2148. y = np.random.randint(0, 16, size=100)
  2149. yf = y.astype(np.float64)
  2150. res = AutoReg(y, 5).fit()
  2151. resf = AutoReg(yf, 5).fit()
  2152. assert_allclose(res.params, resf.params, atol=1e-5)
  2153. assert_allclose(res.bse, resf.bse, atol=1e-5)
  2154. res = ARMA(y, order=(2, 1)).fit(disp=-1)
  2155. resf = ARMA(yf, order=(2, 1)).fit(disp=-1)
  2156. assert_allclose(res.params, resf.params, atol=1e-5)
  2157. assert_allclose(res.bse, resf.bse, atol=1e-5)
  2158. res = ARIMA(y.cumsum(), order=(1, 1, 1)).fit(disp=-1)
  2159. resf = ARIMA(yf.cumsum(), order=(1, 1, 1)).fit(disp=-1)
  2160. assert_allclose(res.params, resf.params, rtol=1e-6, atol=1e-5)
  2161. assert_allclose(res.bse, resf.bse, rtol=1e-6, atol=1e-5)
  2162. def test_multidim_endog(reset_randomstate):
  2163. y = np.random.randn(1000, 2)
  2164. with pytest.raises(ValueError):
  2165. ARMA(y, order=(1, 1))
  2166. y = np.random.randn(1000, 1, 2)
  2167. with pytest.raises(ValueError):
  2168. ARMA(y, order=(1, 1))
  2169. y = np.random.randn(1, 1, 1000)
  2170. with pytest.raises(ValueError):
  2171. ARMA(y, order=(1, 1))
  2172. def test_arima_no_full_output():
  2173. # GH 2752
  2174. endog = y_arma[:, 6]
  2175. mod = ARIMA(endog, (1, 0, 1))
  2176. res = mod.fit(trend="c", disp=-1, full_output=False)
  2177. assert res.mle_retvals is None
  2178. def test_arima_summary_no_lags(reset_randomstate):
  2179. y_train = pd.Series(np.random.randn(1000), name='y_train')
  2180. x_train = pd.DataFrame(np.random.randn(1000, 2), columns=['x1', 'x2'])
  2181. res = ARIMA(endog=y_train.values, exog=x_train,
  2182. order=(0, 1, 0)).fit(disp=-1)
  2183. summ = res.summary().as_text()
  2184. assert 'const ' in summ
  2185. assert 'x1 ' in summ
  2186. assert 'x2 ' in summ
  2187. def test_constant_column_trend():
  2188. # GH#3343 analogous to GH#5258 for AR, when the user passes a constant exog
  2189. # and also passes trend="c", raise instead _make_arma_exog used to
  2190. # silently drop a column and return an inconsistenct k_trend value.
  2191. exog = np.array([0.5, 0.5, 0.5, 0.5, 0.5])
  2192. endog = np.array([-0.011866, 0.003380, 0.015357, 0.004451, -0.020889])
  2193. model = ARIMA(endog=endog, order=(1, 1, 0), exog=exog)
  2194. # Fitting with a constant and constant exog raises because of colinearity
  2195. with pytest.raises(ValueError, match="x contains a constant"):
  2196. model.fit(trend="c", disp=-1)
  2197. # FIXME: calling model.fit(trend="nc") raises for orthogonal reasons
  2198. def test_arima_d0_forecast(reset_randomstate):
  2199. arparams = np.array([1, -.75])
  2200. maparams = np.array([1, .65])
  2201. y = arma_generate_sample(arparams, maparams, 250)
  2202. mod = ARIMA(y, (1, 0, 1))
  2203. assert isinstance(mod, ARMA)
  2204. res = mod.fit(disp=-1)
  2205. pred = res.predict(start=200, end=300)
  2206. pred_typ = res.predict(start=200, end=300, typ='linear')
  2207. assert_allclose(pred, pred_typ)
  2208. pred_mod = mod.predict(res.params, start=200, end=300, typ='linear')
  2209. assert_allclose(pred, pred_mod)
  2210. def test_predict_exog_adjust(reset_randomstate):
  2211. df = pd.DataFrame(np.random.standard_normal((1000, 2)),
  2212. index=pd.date_range('31-12-2000', periods=1000),
  2213. columns=['y', 'x'])
  2214. y = df.y.iloc[:750]
  2215. x = df.x.iloc[:750]
  2216. mod = ARMA(y, (2, 0), exog=x)
  2217. res = mod.fit(disp=-1)
  2218. pred_adj = res.predict(end=df.index[-1], exog=df.x)
  2219. pred = res.predict(end=df.index[-1], exog=df.x.iloc[750:])
  2220. assert_allclose(pred, pred_adj)
  2221. def test_predict_wrong_exog_len(reset_randomstate):
  2222. df = pd.DataFrame(np.random.standard_normal((1000, 2)),
  2223. index=pd.date_range('31-12-2000', periods=1000),
  2224. columns=['y', 'x'])
  2225. y = df.y.iloc[:750]
  2226. x = df.x.iloc[:750]
  2227. mod = ARMA(y, (2, 0), exog=x)
  2228. res = mod.fit(disp=-1)
  2229. with pytest.warns(SpecificationWarning, match='of observations in exog'):
  2230. res.predict(end=df.index[-1], exog=df.x.iloc[500:])
  2231. def test_invalid_method(reset_randomstate):
  2232. ar, ma = [1, -0.5], [1., 0.4]
  2233. x = fa.ArmaFft(ar, ma, 40).generate_sample(nsample=1000, burnin=1000)
  2234. mod = ARMA(x, (1, 1))
  2235. with pytest.raises(ValueError, match='Method unknown not understood'):
  2236. mod.fit(method='unknown')
  2237. @pytest.mark.matplotlib
  2238. def test_summary_plots_no_dates(reset_randomstate, close_figures):
  2239. ar, ma = [1, -0.5], [1., 0.4]
  2240. x = fa.ArmaFft(ar, ma, 40).generate_sample(nsample=1000, burnin=1000)
  2241. res = ARMA(x, (1, 1)).fit(disp=-1)
  2242. res.plot_predict()
  2243. res.plot_predict(plot_insample=True)
  2244. assert isinstance(res.summary().as_text(), str)
  2245. assert isinstance(res.summary2().as_text(), str)
  2246. def test_predict_exog_missing():
  2247. dta = load_macrodata_pandas().data
  2248. res = ARMA(dta['realcons'], (1, 1), dta['m1']).fit(disp=-1)
  2249. with pytest.raises(ValueError, match='You must provide exog for ARMAX'):
  2250. res.predict(start=203, end=210)
  2251. def test_too_few_obs_starting_values(reset_randomstate):
  2252. ar, ma = [1, -0.5], [1., 0.4]
  2253. x = fa.ArmaFft(ar, ma, 40).generate_sample(nsample=5, burnin=100)
  2254. with pytest.raises(ValueError, match='Insufficient degrees of freedom'):
  2255. ARMA(x, (1, 6)).fit(disp=-1)
  2256. def test_repeated_fit_error(reset_randomstate):
  2257. ar, ma = [1, -0.5], [1., 0.4]
  2258. x = fa.ArmaFft(ar, ma, 40).generate_sample(nsample=1000, burnin=100)
  2259. arma = ARMA(x, (1, 1))
  2260. arma.fit(disp=-1)
  2261. with pytest.raises(RuntimeError):
  2262. arma.fit(trend='nc')
  2263. with pytest.raises(RuntimeError):
  2264. arma.fit(method='css')
  2265. def test_predict_no_fit():
  2266. ar, ma = [1, -0.5], [1., 0.4]
  2267. x = fa.ArmaFft(ar, ma, 40).generate_sample(nsample=1000, burnin=100)
  2268. arma = ARMA(x, (1, 1))
  2269. with pytest.raises(RuntimeError, match='Model must be fit'):
  2270. arma.predict([.1, .1])
  2271. arima = ARIMA(x, (1, 1, 1))
  2272. with pytest.raises(RuntimeError, match='Model must be fit'):
  2273. arima.predict([.1, .1])
  2274. def test_arma_repeated_fit():
  2275. ar, ma = [1, -0.5], [1., 0.4]
  2276. x = fa.ArmaFft(ar, ma, 40).generate_sample(nsample=1000, burnin=100)
  2277. arma = ARMA(x, (1, 1))
  2278. res = arma.fit(trend='c', disp=-1)
  2279. repeat = arma.fit(trend='c', disp=-1)
  2280. rtol = 1e-4 if PLATFORM_WIN32 else 1e-7
  2281. assert_allclose(res.params, repeat.params, rtol=rtol)
  2282. assert isinstance(res.summary().as_text(), str)
  2283. assert isinstance(repeat.summary().as_text(), str)
  2284. def test_arima_repeated_fit(reset_randomstate):
  2285. ar, ma = [1, -0.5], [1., 0.4]
  2286. n = 1000
  2287. arma_fft = fa.ArmaFft(ar, ma, 40)
  2288. x = 4 * np.arange(n) + arma_fft.generate_sample(nsample=n, burnin=100)
  2289. arma = ARIMA(x, (1, 1, 1))
  2290. res = arma.fit(trend='c', disp=-1)
  2291. repeat = arma.fit(trend='c', disp=-1)
  2292. tol = 1e-3 if PLATFORM_WIN32 else 1e-6
  2293. assert_allclose(res.params, repeat.params, atol=tol, rtol=tol)
  2294. assert isinstance(res.summary().as_text(), str)
  2295. assert isinstance(repeat.summary().as_text(), str)
  2296. def test_arma_ttest_ftest():
  2297. # GH947
  2298. ar, ma = [1, -0.5], [1., 0.4]
  2299. x = 2 + fa.ArmaFft(ar, ma, 40).generate_sample(nsample=1000, burnin=100)
  2300. arma = ARMA(x, (1, 1))
  2301. res = arma.fit(disp=-1)
  2302. r = np.array([0, 1, 1])
  2303. t_test = res.t_test(r)
  2304. c = res.cov_params()
  2305. v = r @ c @ r
  2306. numerator = res.params @ r
  2307. denom = np.sqrt(v)
  2308. stat = numerator / denom
  2309. assert_allclose(stat, t_test.statistic)
  2310. f_test = res.f_test(r)
  2311. assert_allclose(stat**2, f_test.statistic)
  2312. def test_nan_exog_arima_d1():
  2313. # GH2746, verify raises on missing
  2314. df = pd.DataFrame(np.random.standard_normal((1000, 2)),
  2315. index=pd.date_range('31-12-2000', periods=1000),
  2316. columns=['y', 'x'])
  2317. y = df.y.iloc[:750]
  2318. x = df.x.iloc[:750]
  2319. x.iloc[0] = np.nan
  2320. with pytest.raises(MissingDataError):
  2321. ARIMA(y, (0, 1, 1), exog=x)
  2322. def test_arma_exog_const_trend_nc(reset_randomstate):
  2323. x = np.ones((1000, 1))
  2324. y = arma_generate_sample([1, -1.2, 0.5], [1, 0.7], 1000,
  2325. distrvs=np.random.standard_normal)
  2326. y += 1
  2327. res = ARMA(y, (2, 1), exog=x).fit(trend='nc', disp=-1)
  2328. assert res.params.shape[0] == 4