statsforecast [Model] Arima is getting FloatingPointError

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.<locals>.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__.<locals>.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.

Sep 18 '23 20:09 Timuchin

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

Dec 18 '23 21:12 jmoralez

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.

Jan 18 '24 04:01 github-actions[bot]

@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

Feb 20 '24 07:02 Timuchin

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

Apr 03 '24 19:04 drewbitt

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.

Apr 15 '24 22:04 jmoralez

statsforecast statsforecast copied to clipboard

[Model] Arima is getting FloatingPointError

What happened + What you expected to happen

Versions / Dependencies

Reproduction script

Issue Severity

statsforecast
statsforecast copied to clipboard