PageRenderTime 59ms CodeModel.GetById 20ms RepoModel.GetById 0ms app.codeStats 1ms

/statsmodels/tsa/tests/test_arima.py

https://github.com/danielballan/statsmodels
Python | 1902 lines | 1241 code | 209 blank | 452 comment | 13 complexity | 3a3bac044e13cdeea3436c60cfeede96 MD5 | raw file
Possible License(s): BSD-3-Clause

Large files files are truncated, but you can click here to view the full file

  1. import numpy as np
  2. from numpy.testing import (assert_almost_equal, assert_equal, assert_,
  3. assert_raises, dec)
  4. import statsmodels.sandbox.tsa.fftarma as fa
  5. from statsmodels.tsa.descriptivestats import TsaDescriptive
  6. from statsmodels.tsa.arma_mle import Arma
  7. from statsmodels.tsa.arima_model import ARMA, ARIMA
  8. from statsmodels.tsa.base.datetools import dates_from_range
  9. from results import results_arma, results_arima
  10. import os
  11. from statsmodels.tsa.base import datetools
  12. from statsmodels.tsa.arima_process import arma_generate_sample
  13. import pandas
  14. try:
  15. from statsmodels.tsa.kalmanf import kalman_loglike
  16. fast_kalman = 1
  17. except:
  18. fast_kalman = 0
  19. #NOTE: the KF with complex input returns a different precision for
  20. # the hessian imaginary part, so we use approx_hess and the the
  21. # resulting stats are slightly different.
  22. DECIMAL_4 = 4
  23. DECIMAL_3 = 3
  24. DECIMAL_2 = 2
  25. DECIMAL_1 = 1
  26. current_path = os.path.dirname(os.path.abspath(__file__))
  27. y_arma = np.genfromtxt(open(current_path + '/results/y_arma_data.csv', "rb"),
  28. delimiter=",", skip_header=1, dtype=float)
  29. cpi_dates = dates_from_range('1959Q1', '2009Q3')
  30. cpi_predict_dates = dates_from_range('2009Q3', '2015Q4')
  31. sun_dates = dates_from_range('1700', '2008')
  32. sun_predict_dates = dates_from_range('2008', '2033')
  33. try:
  34. from pandas import DatetimeIndex # pylint: disable-msg=E0611
  35. cpi_dates = DatetimeIndex(cpi_dates, freq='infer')
  36. sun_dates = DatetimeIndex(sun_dates, freq='infer')
  37. cpi_predict_dates = DatetimeIndex(cpi_predict_dates, freq='infer')
  38. sun_predict_dates = DatetimeIndex(sun_predict_dates, freq='infer')
  39. except ImportError:
  40. pass
  41. def test_compare_arma():
  42. #this is a preliminary test to compare arma_kf, arma_cond_ls and arma_cond_mle
  43. #the results returned by the fit methods are incomplete
  44. #for now without random.seed
  45. np.random.seed(9876565)
  46. x = fa.ArmaFft([1, -0.5], [1., 0.4], 40).generate_sample(size=200,
  47. burnin=1000)
  48. # this used kalman filter through descriptive
  49. # d = ARMA(x)
  50. # d.fit((1,1), trend='nc')
  51. # dres = d.res
  52. modkf = ARMA(x, (1,1))
  53. ##rkf = mkf.fit((1,1))
  54. ##rkf.params
  55. reskf = modkf.fit(trend='nc', disp=-1)
  56. dres = reskf
  57. modc = Arma(x)
  58. resls = modc.fit(order=(1,1))
  59. rescm = modc.fit_mle(order=(1,1), start_params=[0.4,0.4, 1.], disp=0)
  60. #decimal 1 corresponds to threshold of 5% difference
  61. #still different sign corrcted
  62. #assert_almost_equal(np.abs(resls[0] / d.params), np.ones(d.params.shape), decimal=1)
  63. assert_almost_equal(resls[0] / dres.params, np.ones(dres.params.shape),
  64. decimal=1)
  65. #rescm also contains variance estimate as last element of params
  66. #assert_almost_equal(np.abs(rescm.params[:-1] / d.params), np.ones(d.params.shape), decimal=1)
  67. assert_almost_equal(rescm.params[:-1] / dres.params, np.ones(dres.params.shape), decimal=1)
  68. #return resls[0], d.params, rescm.params
  69. class CheckArmaResultsMixin(object):
  70. """
  71. res2 are the results from gretl. They are in results/results_arma.
  72. res1 are from statsmodels
  73. """
  74. decimal_params = DECIMAL_4
  75. def test_params(self):
  76. assert_almost_equal(self.res1.params, self.res2.params,
  77. self.decimal_params)
  78. decimal_aic = DECIMAL_4
  79. def test_aic(self):
  80. assert_almost_equal(self.res1.aic, self.res2.aic, self.decimal_aic)
  81. decimal_bic = DECIMAL_4
  82. def test_bic(self):
  83. assert_almost_equal(self.res1.bic, self.res2.bic, self.decimal_bic)
  84. decimal_arroots = DECIMAL_4
  85. def test_arroots(self):
  86. assert_almost_equal(self.res1.arroots, self.res2.arroots,
  87. self.decimal_arroots)
  88. decimal_maroots = DECIMAL_4
  89. def test_maroots(self):
  90. assert_almost_equal(self.res1.maroots, self.res2.maroots,
  91. self.decimal_maroots)
  92. decimal_bse = DECIMAL_2
  93. def test_bse(self):
  94. assert_almost_equal(self.res1.bse, self.res2.bse, self.decimal_bse)
  95. decimal_cov_params = DECIMAL_4
  96. def test_covparams(self):
  97. assert_almost_equal(self.res1.cov_params(), self.res2.cov_params,
  98. self.decimal_cov_params)
  99. decimal_hqic = DECIMAL_4
  100. def test_hqic(self):
  101. assert_almost_equal(self.res1.hqic, self.res2.hqic, self.decimal_hqic)
  102. decimal_llf = DECIMAL_4
  103. def test_llf(self):
  104. assert_almost_equal(self.res1.llf, self.res2.llf, self.decimal_llf)
  105. decimal_resid = DECIMAL_4
  106. def test_resid(self):
  107. assert_almost_equal(self.res1.resid, self.res2.resid,
  108. self.decimal_resid)
  109. decimal_fittedvalues = DECIMAL_4
  110. def test_fittedvalues(self):
  111. assert_almost_equal(self.res1.fittedvalues, self.res2.fittedvalues,
  112. self.decimal_fittedvalues)
  113. decimal_pvalues = DECIMAL_2
  114. def test_pvalues(self):
  115. assert_almost_equal(self.res1.pvalues, self.res2.pvalues,
  116. self.decimal_pvalues)
  117. decimal_t = DECIMAL_2 # only 2 decimal places in gretl output
  118. def test_tvalues(self):
  119. assert_almost_equal(self.res1.tvalues, self.res2.tvalues,
  120. self.decimal_t)
  121. decimal_sigma2 = DECIMAL_4
  122. def test_sigma2(self):
  123. assert_almost_equal(self.res1.sigma2, self.res2.sigma2,
  124. self.decimal_sigma2)
  125. def test_summary(self):
  126. # smoke tests
  127. table = self.res1.summary()
  128. class CheckForecastMixin(object):
  129. decimal_forecast = DECIMAL_4
  130. def test_forecast(self):
  131. assert_almost_equal(self.res1.forecast_res, self.res2.forecast,
  132. self.decimal_forecast)
  133. decimal_forecasterr = DECIMAL_4
  134. def test_forecasterr(self):
  135. assert_almost_equal(self.res1.forecast_err, self.res2.forecasterr,
  136. self.decimal_forecasterr)
  137. class CheckDynamicForecastMixin(object):
  138. decimal_forecast_dyn = 4
  139. def test_dynamic_forecast(self):
  140. assert_almost_equal(self.res1.forecast_res_dyn, self.res2.forecast_dyn,
  141. self.decimal_forecast_dyn)
  142. #def test_forecasterr(self):
  143. # assert_almost_equal(self.res1.forecast_err_dyn,
  144. # self.res2.forecasterr_dyn,
  145. # DECIMAL_4)
  146. class CheckArimaResultsMixin(CheckArmaResultsMixin):
  147. def test_order(self):
  148. assert self.res1.k_diff == self.res2.k_diff
  149. assert self.res1.k_ar == self.res2.k_ar
  150. assert self.res1.k_ma == self.res2.k_ma
  151. decimal_predict_levels = DECIMAL_4
  152. def test_predict_levels(self):
  153. assert_almost_equal(self.res1.predict(typ='levels'), self.res2.linear,
  154. self.decimal_predict_levels)
  155. #NOTE: Ok
  156. class Test_Y_ARMA11_NoConst(CheckArmaResultsMixin, CheckForecastMixin):
  157. @classmethod
  158. def setupClass(cls):
  159. endog = y_arma[:,0]
  160. cls.res1 = ARMA(endog, order=(1,1)).fit(trend='nc', disp=-1)
  161. (cls.res1.forecast_res, cls.res1.forecast_err,
  162. confint) = cls.res1.forecast(10)
  163. cls.res2 = results_arma.Y_arma11()
  164. def test_pickle(self):
  165. from statsmodels.compatnp.py3k import BytesIO
  166. fh = BytesIO()
  167. #test wrapped results load save pickle
  168. self.res1.save(fh)
  169. fh.seek(0,0)
  170. res_unpickled = self.res1.__class__.load(fh)
  171. assert_(type(res_unpickled) is type(self.res1))
  172. #NOTE: Ok
  173. class Test_Y_ARMA14_NoConst(CheckArmaResultsMixin):
  174. @classmethod
  175. def setupClass(cls):
  176. endog = y_arma[:,1]
  177. cls.res1 = ARMA(endog, order=(1,4)).fit(trend='nc', disp=-1)
  178. cls.res2 = results_arma.Y_arma14()
  179. #NOTE: Ok
  180. #can't use class decorators in 2.5....
  181. #@dec.slow
  182. class Test_Y_ARMA41_NoConst(CheckArmaResultsMixin, CheckForecastMixin):
  183. @classmethod
  184. def setupClass(cls):
  185. endog = y_arma[:,2]
  186. cls.res1 = ARMA(endog, order=(4,1)).fit(trend='nc', disp=-1)
  187. (cls.res1.forecast_res, cls.res1.forecast_err,
  188. confint) = cls.res1.forecast(10)
  189. cls.res2 = results_arma.Y_arma41()
  190. cls.decimal_maroots = DECIMAL_3
  191. #NOTE: Ok
  192. class Test_Y_ARMA22_NoConst(CheckArmaResultsMixin):
  193. @classmethod
  194. def setupClass(cls):
  195. endog = y_arma[:,3]
  196. cls.res1 = ARMA(endog, order=(2,2)).fit(trend='nc', disp=-1)
  197. cls.res2 = results_arma.Y_arma22()
  198. #NOTE: Ok
  199. class Test_Y_ARMA50_NoConst(CheckArmaResultsMixin, CheckForecastMixin):
  200. @classmethod
  201. def setupClass(cls):
  202. endog = y_arma[:,4]
  203. cls.res1 = ARMA(endog, order=(5,0)).fit(trend='nc', disp=-1)
  204. (cls.res1.forecast_res, cls.res1.forecast_err,
  205. confint) = cls.res1.forecast(10)
  206. cls.res2 = results_arma.Y_arma50()
  207. #NOTE: Ok
  208. class Test_Y_ARMA02_NoConst(CheckArmaResultsMixin):
  209. @classmethod
  210. def setupClass(cls):
  211. endog = y_arma[:,5]
  212. cls.res1 = ARMA(endog, order=(0,2)).fit(trend='nc', disp=-1)
  213. cls.res2 = results_arma.Y_arma02()
  214. #NOTE: Ok
  215. class Test_Y_ARMA11_Const(CheckArmaResultsMixin, CheckForecastMixin):
  216. @classmethod
  217. def setupClass(cls):
  218. endog = y_arma[:,6]
  219. cls.res1 = ARMA(endog, order=(1,1)).fit(trend="c", disp=-1)
  220. (cls.res1.forecast_res, cls.res1.forecast_err,
  221. confint) = cls.res1.forecast(10)
  222. cls.res2 = results_arma.Y_arma11c()
  223. #NOTE: OK
  224. class Test_Y_ARMA14_Const(CheckArmaResultsMixin):
  225. @classmethod
  226. def setupClass(cls):
  227. endog = y_arma[:,7]
  228. cls.res1 = ARMA(endog, order=(1,4)).fit(trend="c", disp=-1)
  229. cls.res2 = results_arma.Y_arma14c()
  230. #NOTE: Ok
  231. #@dec.slow
  232. class Test_Y_ARMA41_Const(CheckArmaResultsMixin, CheckForecastMixin):
  233. @classmethod
  234. def setupClass(cls):
  235. endog = y_arma[:,8]
  236. cls.res1 = ARMA(endog, order=(4,1)).fit(trend="c", disp=-1)
  237. (cls.res1.forecast_res, cls.res1.forecast_err,
  238. confint) = cls.res1.forecast(10)
  239. cls.res2 = results_arma.Y_arma41c()
  240. cls.decimal_cov_params = DECIMAL_3
  241. cls.decimal_fittedvalues = DECIMAL_3
  242. cls.decimal_resid = DECIMAL_3
  243. cls.decimal_params = DECIMAL_3
  244. if fast_kalman:
  245. cls.decimal_cov_params -= 2
  246. cls.decimal_bse -= 1
  247. #NOTE: Ok
  248. class Test_Y_ARMA22_Const(CheckArmaResultsMixin):
  249. @classmethod
  250. def setupClass(cls):
  251. endog = y_arma[:,9]
  252. cls.res1 = ARMA(endog, order=(2,2)).fit(trend="c", disp=-1)
  253. cls.res2 = results_arma.Y_arma22c()
  254. #NOTE: Ok
  255. class Test_Y_ARMA50_Const(CheckArmaResultsMixin, CheckForecastMixin):
  256. @classmethod
  257. def setupClass(cls):
  258. endog = y_arma[:,10]
  259. cls.res1 = ARMA(endog, order=(5,0)).fit(trend="c", disp=-1)
  260. (cls.res1.forecast_res, cls.res1.forecast_err,
  261. confint) = cls.res1.forecast(10)
  262. cls.res2 = results_arma.Y_arma50c()
  263. #NOTE: Ok
  264. class Test_Y_ARMA02_Const(CheckArmaResultsMixin):
  265. @classmethod
  266. def setupClass(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. #NOTE:
  271. # cov_params and tvalues are off still but not as much vs. R
  272. class Test_Y_ARMA11_NoConst_CSS(CheckArmaResultsMixin):
  273. @classmethod
  274. def setupClass(cls):
  275. endog = y_arma[:,0]
  276. cls.res1 = ARMA(endog, order=(1,1)).fit(method="css", trend='nc',
  277. disp=-1)
  278. cls.res2 = results_arma.Y_arma11("css")
  279. cls.decimal_t = DECIMAL_1
  280. # better vs. R
  281. class Test_Y_ARMA14_NoConst_CSS(CheckArmaResultsMixin):
  282. @classmethod
  283. def setupClass(cls):
  284. endog = y_arma[:,1]
  285. cls.res1 = ARMA(endog, order=(1,4)).fit(method="css", trend='nc',
  286. disp=-1)
  287. cls.res2 = results_arma.Y_arma14("css")
  288. cls.decimal_fittedvalues = DECIMAL_3
  289. cls.decimal_resid = DECIMAL_3
  290. cls.decimal_t = DECIMAL_1
  291. #NOTE: Ok
  292. #NOTE:
  293. # bse, etc. better vs. R
  294. # maroot is off because maparams is off a bit (adjust tolerance?)
  295. class Test_Y_ARMA41_NoConst_CSS(CheckArmaResultsMixin):
  296. @classmethod
  297. def setupClass(cls):
  298. endog = y_arma[:,2]
  299. cls.res1 = ARMA(endog, order=(4,1)).fit(method="css", trend='nc',
  300. disp=-1)
  301. cls.res2 = results_arma.Y_arma41("css")
  302. cls.decimal_t = DECIMAL_1
  303. cls.decimal_pvalues = 0
  304. cls.decimal_cov_params = DECIMAL_3
  305. cls.decimal_maroots = DECIMAL_1
  306. #NOTE: Ok
  307. #same notes as above
  308. class Test_Y_ARMA22_NoConst_CSS(CheckArmaResultsMixin):
  309. @classmethod
  310. def setupClass(cls):
  311. endog = y_arma[:,3]
  312. cls.res1 = ARMA(endog, order=(2,2)).fit(method="css", trend='nc',
  313. disp=-1)
  314. cls.res2 = results_arma.Y_arma22("css")
  315. cls.decimal_t = DECIMAL_1
  316. cls.decimal_resid = DECIMAL_3
  317. cls.decimal_pvalues = DECIMAL_1
  318. cls.decimal_fittedvalues = DECIMAL_3
  319. #NOTE: Ok
  320. #NOTE: gretl just uses least squares for AR CSS
  321. # so BIC, etc. is
  322. # -2*res1.llf + np.log(nobs)*(res1.q+res1.p+res1.k)
  323. # with no adjustment for p and no extra sigma estimate
  324. #NOTE: so our tests use x-12 arima results which agree with us and are
  325. # consistent with the rest of the models
  326. class Test_Y_ARMA50_NoConst_CSS(CheckArmaResultsMixin):
  327. @classmethod
  328. def setupClass(cls):
  329. endog = y_arma[:,4]
  330. cls.res1 = ARMA(endog, order=(5,0)).fit(method="css", trend='nc',
  331. disp=-1)
  332. cls.res2 = results_arma.Y_arma50("css")
  333. cls.decimal_t = 0
  334. cls.decimal_llf = DECIMAL_1 # looks like rounding error?
  335. #NOTE: ok
  336. class Test_Y_ARMA02_NoConst_CSS(CheckArmaResultsMixin):
  337. @classmethod
  338. def setupClass(cls):
  339. endog = y_arma[:,5]
  340. cls.res1 = ARMA(endog, order=(0,2)).fit(method="css", trend='nc',
  341. disp=-1)
  342. cls.res2 = results_arma.Y_arma02("css")
  343. #NOTE: Ok
  344. #NOTE: our results are close to --x-12-arima option and R
  345. class Test_Y_ARMA11_Const_CSS(CheckArmaResultsMixin):
  346. @classmethod
  347. def setupClass(cls):
  348. endog = y_arma[:,6]
  349. cls.res1 = ARMA(endog, order=(1,1)).fit(trend="c", method="css",
  350. disp=-1)
  351. cls.res2 = results_arma.Y_arma11c("css")
  352. cls.decimal_params = DECIMAL_3
  353. cls.decimal_cov_params = DECIMAL_3
  354. cls.decimal_t = DECIMAL_1
  355. #NOTE: Ok
  356. class Test_Y_ARMA14_Const_CSS(CheckArmaResultsMixin):
  357. @classmethod
  358. def setupClass(cls):
  359. endog = y_arma[:,7]
  360. cls.res1 = ARMA(endog, order=(1,4)).fit(trend="c", method="css",
  361. disp=-1)
  362. cls.res2 = results_arma.Y_arma14c("css")
  363. cls.decimal_t = DECIMAL_1
  364. cls.decimal_pvalues = DECIMAL_1
  365. #NOTE: Ok
  366. class Test_Y_ARMA41_Const_CSS(CheckArmaResultsMixin):
  367. @classmethod
  368. def setupClass(cls):
  369. endog = y_arma[:,8]
  370. cls.res1 = ARMA(endog, order=(4,1)).fit(trend="c", method="css",
  371. disp=-1)
  372. cls.res2 = results_arma.Y_arma41c("css")
  373. cls.decimal_t = DECIMAL_1
  374. cls.decimal_cov_params = DECIMAL_1
  375. cls.decimal_maroots = DECIMAL_3
  376. cls.decimal_bse = DECIMAL_1
  377. #NOTE: Ok
  378. class Test_Y_ARMA22_Const_CSS(CheckArmaResultsMixin):
  379. @classmethod
  380. def setupClass(cls):
  381. endog = y_arma[:,9]
  382. cls.res1 = ARMA(endog, order=(2,2)).fit(trend="c", method="css",
  383. disp=-1)
  384. cls.res2 = results_arma.Y_arma22c("css")
  385. cls.decimal_t = 0
  386. cls.decimal_pvalues = DECIMAL_1
  387. #NOTE: Ok
  388. class Test_Y_ARMA50_Const_CSS(CheckArmaResultsMixin):
  389. @classmethod
  390. def setupClass(cls):
  391. endog = y_arma[:,10]
  392. cls.res1 = ARMA(endog, order=(5,0)).fit(trend="c", method="css",
  393. disp=-1)
  394. cls.res2 = results_arma.Y_arma50c("css")
  395. cls.decimal_t = DECIMAL_1
  396. cls.decimal_params = DECIMAL_3
  397. cls.decimal_cov_params = DECIMAL_2
  398. #NOTE: Ok
  399. class Test_Y_ARMA02_Const_CSS(CheckArmaResultsMixin):
  400. @classmethod
  401. def setupClass(cls):
  402. endog = y_arma[:,11]
  403. cls.res1 = ARMA(endog, order=(0,2)).fit(trend="c", method="css",
  404. disp=-1)
  405. cls.res2 = results_arma.Y_arma02c("css")
  406. def test_reset_trend():
  407. endog = y_arma[:,0]
  408. mod = ARMA(endog, order=(1,1))
  409. res1 = mod.fit(trend="c", disp=-1)
  410. res2 = mod.fit(trend="nc", disp=-1)
  411. assert_equal(len(res1.params), len(res2.params)+1)
  412. #@dec.slow
  413. def test_start_params_bug():
  414. data = np.array([1368., 1187, 1090, 1439, 2362, 2783, 2869, 2512, 1804,
  415. 1544, 1028, 869, 1737, 2055, 1947, 1618, 1196, 867, 997, 1862, 2525,
  416. 3250, 4023, 4018, 3585, 3004, 2500, 2441, 2749, 2466, 2157, 1847, 1463,
  417. 1146, 851, 993, 1448, 1719, 1709, 1455, 1950, 1763, 2075, 2343, 3570,
  418. 4690, 3700, 2339, 1679, 1466, 998, 853, 835, 922, 851, 1125, 1299, 1105,
  419. 860, 701, 689, 774, 582, 419, 846, 1132, 902, 1058, 1341, 1551, 1167,
  420. 975, 786, 759, 751, 649, 876, 720, 498, 553, 459, 543, 447, 415, 377,
  421. 373, 324, 320, 306, 259, 220, 342, 558, 825, 994, 1267, 1473, 1601,
  422. 1896, 1890, 2012, 2198, 2393, 2825, 3411, 3406, 2464, 2891, 3685, 3638,
  423. 3746, 3373, 3190, 2681, 2846, 4129, 5054, 5002, 4801, 4934, 4903, 4713,
  424. 4745, 4736, 4622, 4642, 4478, 4510, 4758, 4457, 4356, 4170, 4658, 4546,
  425. 4402, 4183, 3574, 2586, 3326, 3948, 3983, 3997, 4422, 4496, 4276, 3467,
  426. 2753, 2582, 2921, 2768, 2789, 2824, 2482, 2773, 3005, 3641, 3699, 3774,
  427. 3698, 3628, 3180, 3306, 2841, 2014, 1910, 2560, 2980, 3012, 3210, 3457,
  428. 3158, 3344, 3609, 3327, 2913, 2264, 2326, 2596, 2225, 1767, 1190, 792,
  429. 669, 589, 496, 354, 246, 250, 323, 495, 924, 1536, 2081, 2660, 2814, 2992,
  430. 3115, 2962, 2272, 2151, 1889, 1481, 955, 631, 288, 103, 60, 82, 107, 185,
  431. 618, 1526, 2046, 2348, 2584, 2600, 2515, 2345, 2351, 2355, 2409, 2449,
  432. 2645, 2918, 3187, 2888, 2610, 2740, 2526, 2383, 2936, 2968, 2635, 2617,
  433. 2790, 3906, 4018, 4797, 4919, 4942, 4656, 4444, 3898, 3908, 3678, 3605,
  434. 3186, 2139, 2002, 1559, 1235, 1183, 1096, 673, 389, 223, 352, 308, 365,
  435. 525, 779, 894, 901, 1025, 1047, 981, 902, 759, 569, 519, 408, 263, 156,
  436. 72, 49, 31, 41, 192, 423, 492, 552, 564, 723, 921, 1525, 2768, 3531, 3824,
  437. 3835, 4294, 4533, 4173, 4221, 4064, 4641, 4685, 4026, 4323, 4585, 4836,
  438. 4822, 4631, 4614, 4326, 4790, 4736, 4104, 5099, 5154, 5121, 5384, 5274,
  439. 5225, 4899, 5382, 5295, 5349, 4977, 4597, 4069, 3733, 3439, 3052, 2626,
  440. 1939, 1064, 713, 916, 832, 658, 817, 921, 772, 764, 824, 967, 1127, 1153,
  441. 824, 912, 957, 990, 1218, 1684, 2030, 2119, 2233, 2657, 2652, 2682, 2498,
  442. 2429, 2346, 2298, 2129, 1829, 1816, 1225, 1010, 748, 627, 469, 576, 532,
  443. 475, 582, 641, 605, 699, 680, 714, 670, 666, 636, 672, 679, 446, 248, 134,
  444. 160, 178, 286, 413, 676, 1025, 1159, 952, 1398, 1833, 2045, 2072, 1798,
  445. 1799, 1358, 727, 353, 347, 844, 1377, 1829, 2118, 2272, 2745, 4263, 4314,
  446. 4530, 4354, 4645, 4547, 5391, 4855, 4739, 4520, 4573, 4305, 4196, 3773,
  447. 3368, 2596, 2596, 2305, 2756, 3747, 4078, 3415, 2369, 2210, 2316, 2263,
  448. 2672, 3571, 4131, 4167, 4077, 3924, 3738, 3712, 3510, 3182, 3179, 2951,
  449. 2453, 2078, 1999, 2486, 2581, 1891, 1997, 1366, 1294, 1536, 2794, 3211,
  450. 3242, 3406, 3121, 2425, 2016, 1787, 1508, 1304, 1060, 1342, 1589, 2361,
  451. 3452, 2659, 2857, 3255, 3322, 2852, 2964, 3132, 3033, 2931, 2636, 2818,
  452. 3310, 3396, 3179, 3232, 3543, 3759, 3503, 3758, 3658, 3425, 3053, 2620,
  453. 1837, 923, 712, 1054, 1376, 1556, 1498, 1523, 1088, 728, 890, 1413, 2524,
  454. 3295, 4097, 3993, 4116, 3874, 4074, 4142, 3975, 3908, 3907, 3918, 3755,
  455. 3648, 3778, 4293, 4385, 4360, 4352, 4528, 4365, 3846, 4098, 3860, 3230,
  456. 2820, 2916, 3201, 3721, 3397, 3055, 2141, 1623, 1825, 1716, 2232, 2939,
  457. 3735, 4838, 4560, 4307, 4975, 5173, 4859, 5268, 4992, 5100, 5070, 5270,
  458. 4760, 5135, 5059, 4682, 4492, 4933, 4737, 4611, 4634, 4789, 4811, 4379,
  459. 4689, 4284, 4191, 3313, 2770, 2543, 3105, 2967, 2420, 1996, 2247, 2564,
  460. 2726, 3021, 3427, 3509, 3759, 3324, 2988, 2849, 2340, 2443, 2364, 1252,
  461. 623, 742, 867, 684, 488, 348, 241, 187, 279, 355, 423, 678, 1375, 1497,
  462. 1434, 2116, 2411, 1929, 1628, 1635, 1609, 1757, 2090, 2085, 1790, 1846,
  463. 2038, 2360, 2342, 2401, 2920, 3030, 3132, 4385, 5483, 5865, 5595, 5485,
  464. 5727, 5553, 5560, 5233, 5478, 5159, 5155, 5312, 5079, 4510, 4628, 4535,
  465. 3656, 3698, 3443, 3146, 2562, 2304, 2181, 2293, 1950, 1930, 2197, 2796,
  466. 3441, 3649, 3815, 2850, 4005, 5305, 5550, 5641, 4717, 5131, 2831, 3518,
  467. 3354, 3115, 3515, 3552, 3244, 3658, 4407, 4935, 4299, 3166, 3335, 2728,
  468. 2488, 2573, 2002, 1717, 1645, 1977, 2049, 2125, 2376, 2551, 2578, 2629,
  469. 2750, 3150, 3699, 4062, 3959, 3264, 2671, 2205, 2128, 2133, 2095, 1964,
  470. 2006, 2074, 2201, 2506, 2449, 2465, 2064, 1446, 1382, 983, 898, 489, 319,
  471. 383, 332, 276, 224, 144, 101, 232, 429, 597, 750, 908, 960, 1076, 951,
  472. 1062, 1183, 1404, 1391, 1419, 1497, 1267, 963, 682, 777, 906, 1149, 1439,
  473. 1600, 1876, 1885, 1962, 2280, 2711, 2591, 2411])
  474. res = ARMA(data, order=(4,1)).fit(disp=-1)
  475. class Test_ARIMA101(CheckArmaResultsMixin):
  476. # just make sure this works
  477. @classmethod
  478. def setupClass(cls):
  479. endog = y_arma[:,6]
  480. cls.res1 = ARIMA(endog, (1,0,1)).fit(trend="c", disp=-1)
  481. (cls.res1.forecast_res, cls.res1.forecast_err,
  482. confint) = cls.res1.forecast(10)
  483. cls.res2 = results_arma.Y_arma11c()
  484. cls.res2.k_diff = 0
  485. cls.res2.k_ar = 1
  486. cls.res2.k_ma = 1
  487. class Test_ARIMA111(CheckArimaResultsMixin, CheckForecastMixin,
  488. CheckDynamicForecastMixin):
  489. @classmethod
  490. def setupClass(cls):
  491. from statsmodels.datasets.macrodata import load
  492. cpi = load().data['cpi']
  493. cls.res1 = ARIMA(cpi, (1,1,1)).fit(disp=-1)
  494. cls.res2 = results_arima.ARIMA111()
  495. # make sure endog names changes to D.cpi
  496. cls.decimal_llf = 3
  497. cls.decimal_aic = 3
  498. cls.decimal_bic = 3
  499. cls.decimal_cov_params = 2 # this used to be better?
  500. cls.decimal_t = 0
  501. (cls.res1.forecast_res,
  502. cls.res1.forecast_err,
  503. conf_int) = cls.res1.forecast(25)
  504. #cls.res1.forecast_res_dyn = cls.res1.predict(start=164, end=226, typ='levels', dynamic=True)
  505. #TODO: fix the indexing for the end here, I don't think this is right
  506. # if we're going to treat it like indexing
  507. # the forecast from 2005Q1 through 2009Q4 is indices
  508. # 184 through 227 not 226
  509. # note that the first one counts in the count so 164 + 64 is 65
  510. # predictions
  511. cls.res1.forecast_res_dyn = cls.res1.predict(start=164, end=164+63,
  512. typ='levels', dynamic=True)
  513. def test_freq(self):
  514. assert_almost_equal(self.res1.arfreq, [0.0000], 4)
  515. assert_almost_equal(self.res1.mafreq, [0.0000], 4)
  516. class Test_ARIMA111CSS(CheckArimaResultsMixin, CheckForecastMixin,
  517. CheckDynamicForecastMixin):
  518. @classmethod
  519. def setupClass(cls):
  520. from statsmodels.datasets.macrodata import load
  521. cpi = load().data['cpi']
  522. cls.res1 = ARIMA(cpi, (1,1,1)).fit(disp=-1, method='css')
  523. cls.res2 = results_arima.ARIMA111(method='css')
  524. cls.res2.fittedvalues = - cpi[1:-1] + cls.res2.linear
  525. # make sure endog names changes to D.cpi
  526. (cls.res1.forecast_res,
  527. cls.res1.forecast_err,
  528. conf_int) = cls.res1.forecast(25)
  529. cls.decimal_forecast = 2
  530. cls.decimal_forecast_dyn = 2
  531. cls.decimal_forecasterr = 3
  532. cls.res1.forecast_res_dyn = cls.res1.predict(start=164, end=164+63,
  533. typ='levels', dynamic=True)
  534. # precisions
  535. cls.decimal_arroots = 3
  536. cls.decimal_cov_params = 3
  537. cls.decimal_hqic = 3
  538. cls.decimal_maroots = 3
  539. cls.decimal_t = 1
  540. cls.decimal_fittedvalues = 2 # because of rounding when copying
  541. cls.decimal_resid = 2
  542. #cls.decimal_llf = 3
  543. #cls.decimal_aic = 3
  544. #cls.decimal_bic = 3
  545. cls.decimal_predict_levels = DECIMAL_2
  546. class Test_ARIMA112CSS(CheckArimaResultsMixin):
  547. @classmethod
  548. def setupClass(cls):
  549. from statsmodels.datasets.macrodata import load
  550. cpi = load().data['cpi']
  551. cls.res1 = ARIMA(cpi, (1,1,2)).fit(disp=-1, method='css',
  552. start_params = [.905322, -.692425, 1.07366,
  553. 0.172024])
  554. cls.res2 = results_arima.ARIMA112(method='css')
  555. cls.res2.fittedvalues = - cpi[1:-1] + cls.res2.linear
  556. # make sure endog names changes to D.cpi
  557. cls.decimal_llf = 3
  558. cls.decimal_aic = 3
  559. cls.decimal_bic = 3
  560. #(cls.res1.forecast_res,
  561. # cls.res1.forecast_err,
  562. # conf_int) = cls.res1.forecast(25)
  563. #cls.res1.forecast_res_dyn = cls.res1.predict(start=164, end=226, typ='levels', dynamic=True)
  564. #TODO: fix the indexing for the end here, I don't think this is right
  565. # if we're going to treat it like indexing
  566. # the forecast from 2005Q1 through 2009Q4 is indices
  567. # 184 through 227 not 226
  568. # note that the first one counts in the count so 164 + 64 is 65
  569. # predictions
  570. #cls.res1.forecast_res_dyn = self.predict(start=164, end=164+63,
  571. # typ='levels', dynamic=True)
  572. # since we got from gretl don't have linear prediction in differences
  573. cls.decimal_arroots = 3
  574. cls.decimal_maroots = 2
  575. cls.decimal_t = 1
  576. cls.decimal_resid = 2
  577. cls.decimal_fittedvalues = 3
  578. cls.decimal_predict_levels = DECIMAL_3
  579. def test_freq(self):
  580. assert_almost_equal(self.res1.arfreq, [0.5000], 4)
  581. assert_almost_equal(self.res1.mafreq, [0.5000, 0.5000], 4)
  582. #class Test_ARIMADates(CheckArmaResults, CheckForecast, CheckDynamicForecast):
  583. # @classmethod
  584. # def setupClass(cls):
  585. # from statsmodels.datasets.macrodata import load
  586. # from statsmodels.tsa.datetools import dates_from_range
  587. #
  588. # cpi = load().data['cpi']
  589. # dates = dates_from_range('1959q1', length=203)
  590. # cls.res1 = ARIMA(cpi, dates=dates, freq='Q').fit(order=(1,1,1), disp=-1)
  591. # cls.res2 = results_arima.ARIMA111()
  592. # # make sure endog names changes to D.cpi
  593. # cls.decimal_llf = 3
  594. # cls.decimal_aic = 3
  595. # cls.decimal_bic = 3
  596. # (cls.res1.forecast_res,
  597. # cls.res1.forecast_err,
  598. # conf_int) = cls.res1.forecast(25)
  599. def test_arima_predict_mle_dates():
  600. from statsmodels.datasets.macrodata import load
  601. cpi = load().data['cpi']
  602. res1 = ARIMA(cpi, (4,1,1), dates=cpi_dates, freq='Q').fit(disp=-1)
  603. arima_forecasts = np.genfromtxt(open(
  604. current_path + '/results/results_arima_forecasts_all_mle.csv', "rb"),
  605. delimiter=",", skip_header=1, dtype=float)
  606. fc = arima_forecasts[:,0]
  607. fcdyn = arima_forecasts[:,1]
  608. fcdyn2 = arima_forecasts[:,2]
  609. start, end = 2, 51
  610. fv = res1.predict('1959Q3', '1971Q4', typ='levels')
  611. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  612. assert_equal(res1.data.predict_dates, cpi_dates[start:end+1])
  613. start, end = 202, 227
  614. fv = res1.predict('2009Q3', '2015Q4', typ='levels')
  615. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  616. assert_equal(res1.data.predict_dates, cpi_predict_dates)
  617. # make sure dynamic works
  618. start, end = '1960q2', '1971q4'
  619. fv = res1.predict(start, end, dynamic=True, typ='levels')
  620. assert_almost_equal(fv, fcdyn[5:51+1], DECIMAL_4)
  621. start, end = '1965q1', '2015q4'
  622. fv = res1.predict(start, end, dynamic=True, typ='levels')
  623. assert_almost_equal(fv, fcdyn2[24:227+1], DECIMAL_4)
  624. def test_arma_predict_mle_dates():
  625. from statsmodels.datasets.sunspots import load
  626. sunspots = load().data['SUNACTIVITY']
  627. mod = ARMA(sunspots, (9,0), dates=sun_dates, freq='A')
  628. mod.method = 'mle'
  629. assert_raises(ValueError, mod._get_predict_start, *('1701', True))
  630. start, end = 2, 51
  631. _ = mod._get_predict_start('1702', False)
  632. _ = mod._get_predict_end('1751')
  633. assert_equal(mod.data.predict_dates, sun_dates[start:end+1])
  634. start, end = 308, 333
  635. _ = mod._get_predict_start('2008', False)
  636. _ = mod._get_predict_end('2033')
  637. assert_equal(mod.data.predict_dates, sun_predict_dates)
  638. def test_arima_predict_css_dates():
  639. from statsmodels.datasets.macrodata import load
  640. cpi = load().data['cpi']
  641. res1 = ARIMA(cpi, (4,1,1), dates=cpi_dates, freq='Q').fit(disp=-1,
  642. method='css', trend='nc')
  643. params = np.array([ 1.231272508473910,
  644. -0.282516097759915,
  645. 0.170052755782440,
  646. -0.118203728504945,
  647. -0.938783134717947])
  648. arima_forecasts = np.genfromtxt(open(
  649. current_path + '/results/results_arima_forecasts_all_css.csv', "rb"),
  650. delimiter=",", skip_header=1, dtype=float)
  651. fc = arima_forecasts[:,0]
  652. fcdyn = arima_forecasts[:,1]
  653. fcdyn2 = arima_forecasts[:,2]
  654. start, end = 5, 51
  655. fv = res1.model.predict(params, '1960Q2', '1971Q4', typ='levels')
  656. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  657. assert_equal(res1.data.predict_dates, cpi_dates[start:end+1])
  658. start, end = 202, 227
  659. fv = res1.model.predict(params, '2009Q3', '2015Q4', typ='levels')
  660. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  661. assert_equal(res1.data.predict_dates, cpi_predict_dates)
  662. # make sure dynamic works
  663. start, end = 5, 51
  664. fv = res1.model.predict(params, '1960Q2', '1971Q4', typ='levels',
  665. dynamic=True)
  666. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  667. start, end = '1965q1', '2015q4'
  668. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  669. assert_almost_equal(fv, fcdyn2[24:227+1], DECIMAL_4)
  670. def test_arma_predict_css_dates():
  671. from statsmodels.datasets.sunspots import load
  672. sunspots = load().data['SUNACTIVITY']
  673. mod = ARMA(sunspots, (9,0), dates=sun_dates, freq='A')
  674. mod.method = 'css'
  675. assert_raises(ValueError, mod._get_predict_start, *('1701', False))
  676. def test_arima_predict_mle():
  677. from statsmodels.datasets.macrodata import load
  678. cpi = load().data['cpi']
  679. res1 = ARIMA(cpi, (4,1,1)).fit(disp=-1)
  680. # fit the model so that we get correct endog length but use
  681. arima_forecasts = np.genfromtxt(open(
  682. current_path + '/results/results_arima_forecasts_all_mle.csv', "rb"),
  683. delimiter=",", skip_header=1, dtype=float)
  684. fc = arima_forecasts[:,0]
  685. fcdyn = arima_forecasts[:,1]
  686. fcdyn2 = arima_forecasts[:,2]
  687. fcdyn3 = arima_forecasts[:,3]
  688. fcdyn4 = arima_forecasts[:,4]
  689. # 0 indicates the first sample-observation below
  690. # ie., the index after the pre-sample, these are also differenced once
  691. # so the indices are moved back once from the cpi in levels
  692. # start < p, end <p 1959q2 - 1959q4
  693. start, end = 1,3
  694. fv = res1.predict(start, end, typ='levels')
  695. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  696. # start < p, end 0 1959q3 - 1960q1
  697. start, end = 2, 4
  698. fv = res1.predict(start, end, typ='levels')
  699. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  700. # start < p, end >0 1959q3 - 1971q4
  701. start, end = 2, 51
  702. fv = res1.predict(start, end, typ='levels')
  703. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  704. # start < p, end nobs 1959q3 - 2009q3
  705. start, end = 2, 202
  706. fv = res1.predict(start, end, typ='levels')
  707. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  708. # start < p, end >nobs 1959q3 - 2015q4
  709. start, end = 2, 227
  710. fv = res1.predict(start, end, typ='levels')
  711. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  712. # start 0, end >0 1960q1 - 1971q4
  713. start, end = 4, 51
  714. fv = res1.predict(start, end, typ='levels')
  715. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  716. # start 0, end nobs 1960q1 - 2009q3
  717. start, end = 4, 202
  718. fv = res1.predict(start, end, typ='levels')
  719. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  720. # start 0, end >nobs 1960q1 - 2015q4
  721. start, end = 4, 227
  722. fv = res1.predict(start, end, typ='levels')
  723. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  724. # start >p, end >0 1965q1 - 1971q4
  725. start, end = 24, 51
  726. fv = res1.predict(start, end, typ='levels')
  727. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  728. # start >p, end nobs 1965q1 - 2009q3
  729. start, end = 24, 202
  730. fv = res1.predict(start, end, typ='levels')
  731. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  732. # start >p, end >nobs 1965q1 - 2015q4
  733. start, end = 24, 227
  734. fv = res1.predict(start, end, typ='levels')
  735. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  736. # start nobs, end nobs 2009q3 - 2009q3
  737. #NOTE: raises
  738. #start, end = 202, 202
  739. #fv = res1.predict(start, end, typ='levels')
  740. #assert_almost_equal(fv, [])
  741. # start nobs, end >nobs 2009q3 - 2015q4
  742. start, end = 202, 227
  743. fv = res1.predict(start, end, typ='levels')
  744. assert_almost_equal(fv, fc[start:end+1], DECIMAL_3)
  745. # start >nobs, end >nobs 2009q4 - 2015q4
  746. #NOTE: this raises but shouldn't, dynamic forecasts could start
  747. #one period out
  748. start, end = 203, 227
  749. fv = res1.predict(start, end, typ='levels')
  750. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  751. # defaults
  752. start, end = None, None
  753. fv = res1.predict(start, end, typ='levels')
  754. assert_almost_equal(fv, fc[1:203], DECIMAL_4)
  755. #### Dynamic #####
  756. # start < p, end <p 1959q2 - 1959q4
  757. #NOTE: should raise
  758. #start, end = 1,3
  759. #fv = res1.predict(start, end, dynamic=True, typ='levels')
  760. #assert_almost_equal(fv, arima_forecasts[:,15])
  761. # start < p, end 0 1959q3 - 1960q1
  762. #NOTE: below should raise an error
  763. #start, end = 2, 4
  764. #fv = res1.predict(start, end, dynamic=True, typ='levels')
  765. #assert_almost_equal(fv, fcdyn[5:end+1], DECIMAL_4)
  766. # start < p, end >0 1959q3 - 1971q4
  767. #start, end = 2, 51
  768. #fv = res1.predict(start, end, dynamic=True, typ='levels')
  769. #assert_almost_equal(fv, fcdyn[5:end+1], DECIMAL_4)
  770. ## start < p, end nobs 1959q3 - 2009q3
  771. #start, end = 2, 202
  772. #fv = res1.predict(start, end, dynamic=True, typ='levels')
  773. #assert_almost_equal(fv, fcdyn[5:end+1], DECIMAL_4)
  774. ## start < p, end >nobs 1959q3 - 2015q4
  775. #start, end = 2, 227
  776. #fv = res1.predict(start, end, dynamic=True, typ='levels')
  777. #assert_almost_equal(fv, fcdyn[5:end+1], DECIMAL_4)
  778. # start 0, end >0 1960q1 - 1971q4
  779. start, end = 5, 51
  780. fv = res1.predict(start, end, dynamic=True, typ='levels')
  781. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  782. # start 0, end nobs 1960q1 - 2009q3
  783. start, end = 5, 202
  784. fv = res1.predict(start, end, dynamic=True, typ='levels')
  785. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  786. # start 0, end >nobs 1960q1 - 2015q4
  787. start, end = 5, 227
  788. fv = res1.predict(start, end, dynamic=True, typ='levels')
  789. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  790. # start >p, end >0 1965q1 - 1971q4
  791. start, end = 24, 51
  792. fv = res1.predict(start, end, dynamic=True, typ='levels')
  793. assert_almost_equal(fv, fcdyn2[start:end+1], DECIMAL_4)
  794. # start >p, end nobs 1965q1 - 2009q3
  795. start, end = 24, 202
  796. fv = res1.predict(start, end, dynamic=True, typ='levels')
  797. assert_almost_equal(fv, fcdyn2[start:end+1], DECIMAL_4)
  798. # start >p, end >nobs 1965q1 - 2015q4
  799. start, end = 24, 227
  800. fv = res1.predict(start, end, dynamic=True, typ='levels')
  801. assert_almost_equal(fv, fcdyn2[start:end+1], DECIMAL_4)
  802. # start nobs, end nobs 2009q3 - 2009q3
  803. start, end = 202, 202
  804. fv = res1.predict(start, end, dynamic=True, typ='levels')
  805. assert_almost_equal(fv, fcdyn3[start:end+1], DECIMAL_4)
  806. # start nobs, end >nobs 2009q3 - 2015q4
  807. start, end = 202, 227
  808. fv = res1.predict(start, end, dynamic=True, typ='levels')
  809. assert_almost_equal(fv, fcdyn3[start:end+1], DECIMAL_4)
  810. # start >nobs, end >nobs 2009q4 - 2015q4
  811. start, end = 203, 227
  812. fv = res1.predict(start, end, dynamic=True, typ='levels')
  813. assert_almost_equal(fv, fcdyn4[start:end+1], DECIMAL_4)
  814. # defaults
  815. start, end = None, None
  816. fv = res1.predict(start, end, dynamic=True, typ='levels')
  817. assert_almost_equal(fv, fcdyn[5:203], DECIMAL_4)
  818. def _check_start(model, given, expected, dynamic):
  819. start = model._get_predict_start(given, dynamic)
  820. assert_equal(start, expected)
  821. def _check_end(model, given, end_expect, out_of_sample_expect):
  822. end, out_of_sample = model._get_predict_end(given)
  823. assert_equal((end, out_of_sample), (end_expect, out_of_sample_expect))
  824. def test_arma_predict_indices():
  825. from statsmodels.datasets.sunspots import load
  826. sunspots = load().data['SUNACTIVITY']
  827. model = ARMA(sunspots, (9,0), dates=sun_dates, freq='A')
  828. model.method = 'mle'
  829. # raises - pre-sample + dynamic
  830. assert_raises(ValueError, model._get_predict_start, *(0, True))
  831. assert_raises(ValueError, model._get_predict_start, *(8, True))
  832. assert_raises(ValueError, model._get_predict_start, *('1700', True))
  833. assert_raises(ValueError, model._get_predict_start, *('1708', True))
  834. # raises - start out of sample
  835. assert_raises(ValueError, model._get_predict_start, *(311, True))
  836. assert_raises(ValueError, model._get_predict_start, *(311, False))
  837. assert_raises(ValueError, model._get_predict_start, *('2010', True))
  838. assert_raises(ValueError, model._get_predict_start, *('2010', False))
  839. # works - in-sample
  840. # None
  841. # given, expected, dynamic
  842. start_test_cases = [
  843. (None, 9, True),
  844. # all start get moved back by k_diff
  845. (9, 9, True),
  846. (10, 10, True),
  847. # what about end of sample start - last value is first
  848. # forecast
  849. (309, 309, True),
  850. (308, 308, True),
  851. (0, 0, False),
  852. (1, 1, False),
  853. (4, 4, False),
  854. # all start get moved back by k_diff
  855. ('1709', 9, True),
  856. ('1710', 10, True),
  857. # what about end of sample start - last value is first
  858. # forecast
  859. ('2008', 308, True),
  860. ('2009', 309, True),
  861. ('1700', 0, False),
  862. ('1708', 8, False),
  863. ('1709', 9, False),
  864. ]
  865. for case in start_test_cases:
  866. _check_start(*((model,) + case))
  867. # the length of sunspot is 309, so last index is 208
  868. end_test_cases = [(None, 308, 0),
  869. (307, 307, 0),
  870. (308, 308, 0),
  871. (309, 308, 1),
  872. (312, 308, 4),
  873. (51, 51, 0),
  874. (333, 308, 25),
  875. ('2007', 307, 0),
  876. ('2008', 308, 0),
  877. ('2009', 308, 1),
  878. ('2012', 308, 4),
  879. ('1815', 115, 0),
  880. ('2033', 308, 25),
  881. ]
  882. for case in end_test_cases:
  883. _check_end(*((model,)+case))
  884. def test_arima_predict_indices():
  885. from statsmodels.datasets.macrodata import load
  886. cpi = load().data['cpi']
  887. model = ARIMA(cpi, (4,1,1), dates=cpi_dates, freq='Q')
  888. model.method = 'mle'
  889. # starting indices
  890. # raises - pre-sample + dynamic
  891. assert_raises(ValueError, model._get_predict_start, *(0, True))
  892. assert_raises(ValueError, model._get_predict_start, *(4, True))
  893. assert_raises(ValueError, model._get_predict_start, *('1959Q1', True))
  894. assert_raises(ValueError, model._get_predict_start, *('1960Q1', True))
  895. # raises - index differenced away
  896. assert_raises(ValueError, model._get_predict_start, *(0, False))
  897. assert_raises(ValueError, model._get_predict_start, *('1959Q1', False))
  898. # raises - start out of sample
  899. assert_raises(ValueError, model._get_predict_start, *(204, True))
  900. assert_raises(ValueError, model._get_predict_start, *(204, False))
  901. assert_raises(ValueError, model._get_predict_start, *('2010Q1', True))
  902. assert_raises(ValueError, model._get_predict_start, *('2010Q1', False))
  903. # works - in-sample
  904. # None
  905. # given, expected, dynamic
  906. start_test_cases = [
  907. (None, 4, True),
  908. # all start get moved back by k_diff
  909. (5, 4, True),
  910. (6, 5, True),
  911. # what about end of sample start - last value is first
  912. # forecast
  913. (203, 202, True),
  914. (1, 0, False),
  915. (4, 3, False),
  916. (5, 4, False),
  917. # all start get moved back by k_diff
  918. ('1960Q2', 4, True),
  919. ('1960Q3', 5, True),
  920. # what about end of sample start - last value is first
  921. # forecast
  922. ('2009Q4', 202, True),
  923. ('1959Q2', 0, False),
  924. ('1960Q1', 3, False),
  925. ('1960Q2', 4, False),
  926. ]
  927. for case in start_test_cases:
  928. _check_start(*((model,) + case))
  929. # check raises
  930. #TODO: make sure dates are passing through unmolested
  931. #assert_raises(ValueError, model._get_predict_end, ("2001-1-1",))
  932. # the length of diff(cpi) is 202, so last index is 201
  933. end_test_cases = [(None, 201, 0),
  934. (201, 200, 0),
  935. (202, 201, 0),
  936. (203, 201, 1),
  937. (204, 201, 2),
  938. (51, 50, 0),
  939. (164+63, 201, 25),
  940. ('2009Q2', 200, 0),
  941. ('2009Q3', 201, 0),
  942. ('2009Q4', 201, 1),
  943. ('2010Q1', 201, 2),
  944. ('1971Q4', 50, 0),
  945. ('2015Q4', 201, 25),
  946. ]
  947. for case in end_test_cases:
  948. _check_end(*((model,)+case))
  949. # check higher k_diff
  950. model.k_diff = 2
  951. # raises - pre-sample + dynamic
  952. assert_raises(ValueError, model._get_predict_start, *(0, True))
  953. assert_raises(ValueError, model._get_predict_start, *(5, True))
  954. assert_raises(ValueError, model._get_predict_start, *('1959Q1', True))
  955. assert_raises(ValueError, model._get_predict_start, *('1960Q1', True))
  956. # raises - index differenced away
  957. assert_raises(ValueError, model._get_predict_start, *(1, False))
  958. assert_raises(ValueError, model._get_predict_start, *('1959Q2', False))
  959. start_test_cases = [(None, 4, True),
  960. # all start get moved back by k_diff
  961. (6, 4, True),
  962. # what about end of sample start - last value is first
  963. # forecast
  964. (203, 201, True),
  965. (2, 0, False),
  966. (4, 2, False),
  967. (5, 3, False),
  968. ('1960Q3', 4, True),
  969. # what about end of sample start - last value is first
  970. # forecast
  971. ('2009Q4', 201, True),
  972. ('2009Q4', 201, True),
  973. ('1959Q3', 0, False),
  974. ('1960Q1', 2, False),
  975. ('1960Q2', 3, False),
  976. ]
  977. for case in start_test_cases:
  978. _check_start(*((model,)+case))
  979. end_test_cases = [(None, 200, 0),
  980. (201, 199, 0),
  981. (202, 200, 0),
  982. (203, 200, 1),
  983. (204, 200, 2),
  984. (51, 49, 0),
  985. (164+63, 200, 25),
  986. ('2009Q2', 199, 0),
  987. ('2009Q3', 200, 0),
  988. ('2009Q4', 200, 1),
  989. ('2010Q1', 200, 2),
  990. ('1971Q4', 49, 0),
  991. ('2015Q4', 200, 25),
  992. ]
  993. for case in end_test_cases:
  994. _check_end(*((model,)+case))
  995. def test_arima_predict_indices_css():
  996. from statsmodels.datasets.macrodata import load
  997. cpi = load().data['cpi']
  998. #NOTE: Doing no-constant for now to kick the conditional exogenous
  999. #issue 274 down the road
  1000. # go ahead and git the model to set up necessary variables
  1001. model = ARIMA(cpi, (4,1,1))
  1002. model.method = 'css'
  1003. assert_raises(ValueError, model._get_predict_start, *(0, False))
  1004. assert_raises(ValueError, model._get_predict_start, *(0, True))
  1005. assert_raises(ValueError, model._get_predict_start, *(2, False))
  1006. assert_raises(ValueError, model._get_predict_start, *(2, True))
  1007. def test_arima_predict_css():
  1008. from statsmodels.datasets.macrodata import load
  1009. cpi = load().data['cpi']
  1010. #NOTE: Doing no-constant for now to kick the conditional exogenous
  1011. #issue 274 down the road
  1012. # go ahead and git the model to set up necessary variables
  1013. res1 = ARIMA(cpi, (4,1,1)).fit(disp=-1, method="css",
  1014. trend="nc")
  1015. # but use gretl parameters to predict to avoid precision problems
  1016. params = np.array([ 1.231272508473910,
  1017. -0.282516097759915,
  1018. 0.170052755782440,
  1019. -0.118203728504945,
  1020. -0.938783134717947])
  1021. arima_forecasts = np.genfromtxt(open(
  1022. current_path + '/results/results_arima_forecasts_all_css.csv', "rb"),
  1023. delimiter=",", skip_header=1, dtype=float)
  1024. fc = arima_forecasts[:,0]
  1025. fcdyn = arima_forecasts[:,1]
  1026. fcdyn2 = arima_forecasts[:,2]
  1027. fcdyn3 = arima_forecasts[:,3]
  1028. fcdyn4 = arima_forecasts[:,4]
  1029. #NOTE: should raise
  1030. #start, end = 1,3
  1031. #fv = res1.model.predict(params, start, end)
  1032. ## start < p, end 0 1959q3 - 1960q1
  1033. #start, end = 2, 4
  1034. #fv = res1.model.predict(params, start, end)
  1035. ## start < p, end >0 1959q3 - 1971q4
  1036. #start, end = 2, 51
  1037. #fv = res1.model.predict(params, start, end)
  1038. ## start < p, end nobs 1959q3 - 2009q3
  1039. #start, end = 2, 202
  1040. #fv = res1.model.predict(params, start, end)
  1041. ## start < p, end >nobs 1959q3 - 2015q4
  1042. #start, end = 2, 227
  1043. #fv = res1.model.predict(params, start, end)
  1044. # start 0, end >0 1960q1 - 1971q4
  1045. start, end = 5, 51
  1046. fv = res1.model.predict(params, start, end, typ='levels')
  1047. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1048. # start 0, end nobs 1960q1 - 2009q3
  1049. start, end = 5, 202
  1050. fv = res1.model.predict(params, start, end, typ='levels')
  1051. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1052. # start 0, end >nobs 1960q1 - 2015q4
  1053. #TODO: why detoriating precision?
  1054. fv = res1.model.predict(params, start, end, typ='levels')
  1055. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1056. # start >p, end >0 1965q1 - 1971q4
  1057. start, end = 24, 51
  1058. fv = res1.model.predict(params, start, end, typ='levels')
  1059. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1060. # start >p, end nobs 1965q1 - 2009q3
  1061. start, end = 24, 202
  1062. fv = res1.model.predict(params, start, end, typ='levels')
  1063. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1064. # start >p, end >nobs 1965q1 - 2015q4
  1065. start, end = 24, 227
  1066. fv = res1.model.predict(params, start, end, typ='levels')
  1067. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1068. # start nobs, end nobs 2009q3 - 2009q3
  1069. start, end = 202, 202
  1070. fv = res1.model.predict(params, start, end, typ='levels')
  1071. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1072. # start nobs, end >nobs 2009q3 - 2015q4
  1073. start, end = 202, 227
  1074. fv = res1.model.predict(params, start, end, typ='levels')
  1075. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1076. # start >nobs, end >nobs 2009q4 - 2015q4
  1077. start, end = 203, 227
  1078. fv = res1.model.predict(params, start, end, typ='levels')
  1079. assert_almost_equal(fv, fc[start:end+1], DECIMAL_4)
  1080. # defaults
  1081. start, end = None, None
  1082. fv = res1.model.predict(params, start, end, typ='levels')
  1083. assert_almost_equal(fv, fc[5:203], DECIMAL_4)
  1084. #### Dynamic #####
  1085. #NOTE: should raise
  1086. # start < p, end <p 1959q2 - 1959q4
  1087. #start, end = 1,3
  1088. #fv = res1.predict(start, end, dynamic=True)
  1089. # start < p, end 0 1959q3 - 1960q1
  1090. #start, end = 2, 4
  1091. #fv = res1.predict(start, end, dynamic=True)
  1092. ## start < p, end >0 1959q3 - 1971q4
  1093. #start, end = 2, 51
  1094. #fv = res1.predict(start, end, dynamic=True)
  1095. ## start < p, end nobs 1959q3 - 2009q3
  1096. #start, end = 2, 202
  1097. #fv = res1.predict(start, end, dynamic=True)
  1098. ## start < p, end >nobs 1959q3 - 2015q4
  1099. #start, end = 2, 227
  1100. #fv = res1.predict(start, end, dynamic=True)
  1101. # start 0, end >0 1960q1 - 1971q4
  1102. start, end = 5, 51
  1103. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1104. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  1105. # start 0, end nobs 1960q1 - 2009q3
  1106. start, end = 5, 202
  1107. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1108. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  1109. # start 0, end >nobs 1960q1 - 2015q4
  1110. start, end = 5, 227
  1111. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1112. assert_almost_equal(fv, fcdyn[start:end+1], DECIMAL_4)
  1113. # start >p, end >0 1965q1 - 1971q4
  1114. start, end = 24, 51
  1115. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1116. assert_almost_equal(fv, fcdyn2[start:end+1], DECIMAL_4)
  1117. # start >p, end nobs 1965q1 - 2009q3
  1118. start, end = 24, 202
  1119. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1120. assert_almost_equal(fv, fcdyn2[start:end+1], DECIMAL_4)
  1121. # start >p, end >nobs 1965q1 - 2015q4
  1122. start, end = 24, 227
  1123. fv = res1.model.predict(params, start, end, dynamic=True, typ='levels')
  1124. assert_almost_equal(fv

Large files files are truncated, but you can click here to view the full file