[Model] Arima is getting FloatingPointError

Timuchin commented 9 months ago

What happened + What you expected to happen

Hello!

First of all, thank you for your work!

I am comparing ARIMA and other models. Found strange behavior.

MA(1) cannot be fitted and receiving floating point error. As I understand it happens during scipy.minimize.

Interesting enough that ARIMA from statsmodels does not have such problem with this data. Though, it throws warning about Non-invertible starting MA parameters found.

Stacktrace

``` --------------------------------------------------------------------------- FloatingPointError Traceback (most recent call last) [c:\git\forecast-tool\examples\active_listers_example\active_listers.ipynb](file:///C:/git/forecast-tool/examples/active_listers_example/active_listers.ipynb) Cell 200 line 1 [16](vscode-notebook-cell:/c%3A/git/forecast-tool/examples/active_listers_example/active_listers.ipynb#Z1021sZmlsZQ%3D%3D?line=15) eps.append(LOC + 0.99 * eps[i-1] + np.random.normal(loc = np.log(LOC), scale = np.sqrt(SCALE), size = 1)[0]) [17](vscode-notebook-cell:/c%3A/git/forecast-tool/examples/active_listers_example/active_listers.ipynb#Z1021sZmlsZQ%3D%3D?line=16) eps = np.array(eps) ---> [19](vscode-notebook-cell:/c%3A/git/forecast-tool/examples/active_listers_example/active_listers.ipynb#Z1021sZmlsZQ%3D%3D?line=18) arima(x = eps, order = (0, 0, 1)) File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\statsforecast\arima.py:874](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/statsforecast/arima.py:874), in arima(x, order, seasonal, xreg, include_mean, transform_pars, fixed, init, method, SSinit, optim_method, kappa, tol, optim_control) 870 res = OptimResult( 871 True, 0, np.array([]), arma_css_op(np.array([]), x), np.array([]) 872 ) 873 else: --> 874 res = minimize( 875 arma_css_op, 876 init0, 877 args=(x,), 878 method=optim_method, 879 tol=tol, 880 options=optim_control, 881 ) 883 if res.status > 0: 884 warnings.warn( 885 f"possible convergence problem: minimize gave code {res.status}]" 886 ) File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_minimize.py:705](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_minimize.py:705), in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options) 703 res = _minimize_cg(fun, x0, args, jac, callback, **options) 704 elif meth == 'bfgs': --> 705 res = _minimize_bfgs(fun, x0, args, jac, callback, **options) 706 elif meth == 'newton-cg': 707 res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback, 708 **options) File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_optimize.py:1444](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_optimize.py:1444), in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, finite_diff_rel_step, xrtol, **unknown_options) 1441 pk = -np.dot(Hk, gfk) 1442 try: 1443 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \ -> 1444 _line_search_wolfe12(f, myfprime, xk, pk, gfk, 1445 old_fval, old_old_fval, amin=1e-100, amax=1e100) 1446 except _LineSearchError: 1447 # Line search failed to find a better solution. 1448 warnflag = 2 File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_optimize.py:1215](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_optimize.py:1215), in _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval, **kwargs) 1201 """ 1202 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if 1203 suitable step length is not found, and raise an exception if a (...) 1210 1211 """ 1213 extra_condition = kwargs.pop('extra_condition', None) -> 1215 ret = line_search_wolfe1(f, fprime, xk, pk, gfk, 1216 old_fval, old_old_fval, 1217 **kwargs) 1219 if ret[0] is not None and extra_condition is not None: 1220 xp1 = xk + ret[0] * pk File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_linesearch.py:84](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_linesearch.py:84), in line_search_wolfe1(f, fprime, xk, pk, gfk, old_fval, old_old_fval, args, c1, c2, amax, amin, xtol) 80 return np.dot(gval[0], pk) 82 derphi0 = np.dot(gfk, pk) ---> 84 stp, fval, old_fval = scalar_search_wolfe1( 85 phi, derphi, old_fval, old_old_fval, derphi0, 86 c1=c1, c2=c2, amax=amax, amin=amin, xtol=xtol) 88 return stp, fc[0], gc[0], fval, old_fval, gval[0] File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_linesearch.py:161](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_linesearch.py:161), in scalar_search_wolfe1(phi, derphi, phi0, old_phi0, derphi0, c1, c2, amax, amin, xtol) 159 alpha1 = stp 160 phi1 = phi(stp) --> 161 derphi1 = derphi(stp) 162 else: 163 break File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_linesearch.py:78](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_linesearch.py:78), in line_search_wolfe1..derphi(s) 77 def derphi(s): ---> 78 gval[0] = fprime(xk + s*pk, *args) 79 gc[0] += 1 80 return np.dot(gval[0], pk) File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_differentiable_functions.py:273](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_differentiable_functions.py:273), in ScalarFunction.grad(self, x) 271 if not np.array_equal(x, self.x): 272 self._update_x_impl(x) --> 273 self._update_grad() 274 return self.g File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_differentiable_functions.py:256](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_differentiable_functions.py:256), in ScalarFunction._update_grad(self) 254 def _update_grad(self): 255 if not self.g_updated: --> 256 self._update_grad_impl() 257 self.g_updated = True File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_differentiable_functions.py:173](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_differentiable_functions.py:173), in ScalarFunction.__init__..update_grad() 171 self._update_fun() 172 self.ngev += 1 --> 173 self.g = approx_derivative(fun_wrapped, self.x, f0=self.f, 174 **finite_diff_options) File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_numdiff.py:505](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_numdiff.py:505), in approx_derivative(fun, x0, method, rel_step, abs_step, f0, bounds, sparsity, as_linear_operator, args, kwargs) 502 use_one_sided = False 504 if sparsity is None: --> 505 return _dense_difference(fun_wrapped, x0, f0, h, 506 use_one_sided, method) 507 else: 508 if not issparse(sparsity) and len(sparsity) == 2: File [c:\Users\tazekokh\AppData\Local\anaconda3\envs\forecast\Lib\site-packages\scipy\optimize\_numdiff.py:576](file:///C:/Users/tazekokh/AppData/Local/anaconda3/envs/forecast/Lib/site-packages/scipy/optimize/_numdiff.py:576), in _dense_difference(fun, x0, f0, h, use_one_sided, method) 574 x = x0 + h_vecs[i] 575 dx = x[i] - x0[i] # Recompute dx as exactly representable number. --> 576 df = fun(x) - f0 577 elif method == '3-point' and use_one_sided[i]: 578 x1 = x0 + h_vecs[i] FloatingPointError: invalid value encountered in subtract ```

Versions / Dependencies

Dependencies

Windows 10 python=3.11.5 numpy=1.24.4 statsforecast=1.6.0 All conda env dependencies. - adagio=0.2.4=pyhd8ed1ab_0 - antlr4-python3-runtime=4.11.1=py311haa95532_0 - appdirs=1.4.4=pyh9f0ad1d_0 - aws-c-auth=0.7.3=hd125877_3 - aws-c-cal=0.6.2=hfb91821_0 - aws-c-common=0.9.0=hcfcfb64_0 - aws-c-compression=0.2.17=h04c9df6_2 - aws-c-event-stream=0.3.2=h495bb32_0 - aws-c-http=0.7.12=h0890e15_1 - aws-c-io=0.13.32=h83b3346_3 - aws-c-mqtt=0.9.5=h0fd1aac_1 - aws-c-s3=0.3.17=h9f49523_0 - aws-c-sdkutils=0.1.12=h04c9df6_1 - aws-checksums=0.1.17=h04c9df6_1 - aws-crt-cpp=0.23.1=hfe9bf68_1 - aws-sdk-cpp=1.11.156=h03febf0_2 - brotli=1.1.0=hcfcfb64_0 - brotli-bin=1.1.0=hcfcfb64_0 - bzip2=1.0.8=h8ffe710_4 - c-ares=1.19.1=hcfcfb64_0 - ca-certificates=2023.7.22=h56e8100_0 - certifi=2023.7.22=pyhd8ed1ab_0 - colorama=0.4.6=pyhd8ed1ab_0 - contourpy=1.1.1=py311h005e61a_0 - cycler=0.11.0=pyhd8ed1ab_0 - fonttools=4.42.1=py311ha68e1ae_0 - freetype=2.12.1=hdaf720e_2 - fs=2.4.16=pyhd8ed1ab_0 - fsspec=2023.9.1=pyh1a96a4e_0 - fugue=0.8.6=pyhd8ed1ab_0 - fugue-sql-antlr=0.1.7=pyhd8ed1ab_0 - importlib-metadata=6.8.0=pyha770c72_0 - importlib_metadata=6.8.0=hd8ed1ab_0 - intel-openmp=2023.2.0=h57928b3_49496 - jinja2=3.1.2=pyhd8ed1ab_1 - kiwisolver=1.4.5=py311h005e61a_0 - krb5=1.21.2=heb0366b_0 - lcms2=2.15=he9d350c_2 - lerc=4.0.0=h63175ca_0 - libabseil=20230802.0=cxx17_h63175ca_3 - libarrow=13.0.0=h1e3473c_4_cpu - libblas=3.9.0=18_win64_mkl - libbrotlicommon=1.1.0=hcfcfb64_0 - libbrotlidec=1.1.0=hcfcfb64_0 - libbrotlienc=1.1.0=hcfcfb64_0 - libcblas=3.9.0=18_win64_mkl - libcrc32c=1.1.2=h0e60522_0 - libcurl=8.3.0=hd5e4a3a_0 - libdeflate=1.19=hcfcfb64_0 - libevent=2.1.12=h3671451_1 - libexpat=2.5.0=h63175ca_1 - libffi=3.4.2=h8ffe710_5 - libgoogle-cloud=2.12.0=h0a0a397_2 - libgrpc=1.57.0=h550f6bd_1 - libhwloc=2.9.2=default_haede6df_1009 - libiconv=1.17=h8ffe710_0 - libjpeg-turbo=2.1.5.1=hcfcfb64_1 - liblapack=3.9.0=18_win64_mkl - libpng=1.6.39=h19919ed_0 - libprotobuf=4.23.4=hb8276f3_6 - libsqlite=3.43.0=hcfcfb64_0 - libssh2=1.11.0=h7dfc565_0 - libthrift=0.19.0=h06f6336_0 - libtiff=4.6.0=h4554b19_1 - libutf8proc=2.8.0=h82a8f57_0 - libwebp-base=1.3.2=hcfcfb64_0 - libxcb=1.15=hcd874cb_0 - libxml2=2.11.5=hc3477c8_1 - libzlib=1.2.13=hcfcfb64_5 - llvmlite=0.40.1=py311h5bc0dda_0 - lz4-c=1.9.4=hcfcfb64_0 - m2w64-gcc-libgfortran=5.3.0=6 - m2w64-gcc-libs=5.3.0=7 - m2w64-gcc-libs-core=5.3.0=7 - m2w64-gmp=6.1.0=2 - m2w64-libwinpthread-git=5.0.0.4634.697f757=2 - markupsafe=2.1.3=py311ha68e1ae_0 - matplotlib-base=3.7.1=py311h6e989c2_0 - mkl=2022.1.0=h6a75c08_874 - msys2-conda-epoch=20160418=1 - munkres=1.1.4=pyh9f0ad1d_0 - numba=0.57.1=py311h2c0921f_0 - numpy=1.24.4=py311h0b4df5a_0 - openjpeg=2.5.0=h3d672ee_3 - openssl=3.1.2=hcfcfb64_0 - orc=1.9.0=h8dbeef6_2 - packaging=23.1=pyhd8ed1ab_0 - pandas=2.1.0=py311hf63dbb6_0 - patsy=0.5.3=pyhd8ed1ab_0 - pillow=10.0.1=py311hd926f49_0 - pip=23.2.1=pyhd8ed1ab_0 - polars=0.19.3=py311h9dab8b3_0 - pthread-stubs=0.4=hcd874cb_1001 - pthreads-win32=2.9.1=hfa6e2cd_3 - pyarrow=13.0.0=py311h6a6099b_4_cpu - pyparsing=3.1.1=pyhd8ed1ab_0 - python=3.11.5=h2628c8c_0_cpython - python-dateutil=2.8.2=pyhd8ed1ab_0 - python-tzdata=2023.3=pyhd8ed1ab_0 - python_abi=3.11=3_cp311 - pytz=2023.3.post1=pyhd8ed1ab_0 - qpd=0.4.4=pyhd8ed1ab_1 - re2=2023.03.02=hd4eee63_0 - scipy=1.11.2=py311h37ff6ca_1 - setuptools=68.2.2=pyhd8ed1ab_0 - six=1.16.0=pyh6c4a22f_0 - snappy=1.1.10=hfb803bf_0 - sqlglot=18.5.1=pyhd8ed1ab_0 - statsforecast=1.6.0=pyhd8ed1ab_0 - statsmodels=0.14.0=py311h59ca53f_1 - tbb=2021.10.0=h91493d7_0 - tk=8.6.12=h8ffe710_0 - tqdm=4.66.1=pyhd8ed1ab_0 - triad=0.9.1=pyhd8ed1ab_0 - tzdata=2023c=h71feb2d_0 - ucrt=10.0.22621.0=h57928b3_0 - vc=14.3=h64f974e_17 - vc14_runtime=14.36.32532=hdcecf7f_17 - vs2015_runtime=14.36.32532=h05e6639_17 - wheel=0.41.2=pyhd8ed1ab_0 - xorg-libxau=1.0.11=hcd874cb_0 - xorg-libxdmcp=1.1.3=hcd874cb_0 - xz=5.2.6=h8d14728_0 - zipp=3.16.2=pyhd8ed1ab_0 - zstd=1.5.5=h12be248_0

Reproduction script

from statsforecast.arima import arima
import numpy as np

np.random.seed(1234)

SIZE = 1_000
LOC = 1
SCALE = 1

eps = []

for i in range(SIZE):
    if i == 0:
        eps.append(np.random.normal(loc = LOC, scale = SCALE, size = 1)[0])
    elif (i > 0) & (i <= SIZE):
        eps.append(LOC + 0.99 * eps[i-1] + np.random.normal(loc = np.log(LOC), scale = np.sqrt(SCALE), size = 1)[0])

eps = np.array(eps)

arima(x = eps, order = (0, 0, 1))

Issue Severity

Medium: It is a significant difficulty but I can work around it.

jmoralez commented 6 months ago

Hey @Timuchin, thanks for the excellent report. Are you still getting this error? I get it as a warning, i.e. RuntimeWarning: invalid value encountered in subtract df = fun(x) - f0

github-actions[bot] commented 5 months ago

This issue has been automatically closed because it has been awaiting a response for too long. When you have time to to work with the maintainers to resolve this issue, please post a new comment and it will be re-opened. If the issue has been locked for editing by the time you return to it, please open a new issue and reference this one.

Timuchin commented 4 months ago

@jmoralez Hello! After updating statsforecast to 1.7.3 it stopped getting error. But yeah, it throws a warning. I guess the main reason that it cannot converge to stable solution. RuntimeWarning: invalid value encountered in subtract df = fun(x) - f0 UserWarning: possible convergence problem: minimize gave code 2]

It does not seem as a big problem, though quite interesting that MA(1) from statsmodels does not get such error. I guess it successfully converges. As a result coefficient for MA(1) is different. And, what's more important, it calculates AIC criteria. This can can be important further for AutoARIMA.

param	statsforecast	statsmodels tsa	statsmodels statespace
MA	0.911	0.975	0.975
intercept	91.29	91.71	91.71
AIC	nan	7561.34	7561.34

from statsforecast.arima import arima
import numpy as np
import statsmodels.api as sm

np.random.seed(1234)

SIZE = 1_000
LOC = 1
SCALE = 1

eps = []

for i in range(SIZE):
    if i == 0:
        eps.append(np.random.normal(loc = LOC, scale = SCALE, size = 1)[0])
    elif (i > 0) & (i <= SIZE):
        eps.append(LOC + 0.99 * eps[i-1] + np.random.normal(loc = np.log(LOC), scale = np.sqrt(SCALE), size = 1)[0])

eps = np.array(eps)

ma_nixtla = arima(x = eps, order = (0, 0, 1))
ma_stats = sm.tsa.ARIMA(eps, order = (0, 0, 1)).fit()
ma_stats_ss  = sm.tsa.statespace.SARIMAX(eps, order = (0, 0, 1), trend='c').fit()

print(f'''
      MA coef are statsforecast:{ma_nixtla.get('coef').get('ma1')}, statsmodels tsa: {ma_stats.params[1]}, statsmodels statespace: {ma_stats_ss.params[1]}, \n
      Intercept coef are statsforecast:{ma_nixtla.get('coef').get('intercept')}, statsmodels tsa: {ma_stats.params[0]}, statsmodels statespace: {ma_stats_ss.params[0]}, \n
      AIC criteria are statsforecast:{ma_nixtla.get('aic')}, statsmodels tsa: {ma_stats.aic}, statsmodels statespace: {ma_stats_ss.aic}
      ''')

Dependencies

statsmodels==0.14.1

statsforecast==1.7.3

numpy==1.26.3

drewbitt commented 2 months ago

I am also seeing a ton of RuntimeWarning: invalid value encountered in subtract at df = fun(x) - f0

jmoralez commented 2 months ago

Hey @Timuchin. That seems to be due to the default method ('CSS') of the arima function, by setting method='CSS-ML' (the default of statsforecast.models.ARIMA) I get the same result as in statsmodels.

The warnings seem to be coming from the suggestions of the optimization algorithm, which sometimes suggests values that produce Inf in the objective function. I think we can make that return a big float instead.

Nixtla / statsforecast