diff --git a/.gitignore b/.gitignore index 4032fd610..cd5df2718 100644 --- a/.gitignore +++ b/.gitignore @@ -25,3 +25,4 @@ nbs/.last_checked mlruns/ .luarc.json *.so +*.dll diff --git a/experiments/m3/src/evaluation.py b/experiments/m3/src/evaluation.py index 6df6f4350..4a082ecfe 100644 --- a/experiments/m3/src/evaluation.py +++ b/experiments/m3/src/evaluation.py @@ -50,11 +50,11 @@ def main(test: bool = False): time = evaluation.query('metric=="time"').T if test: expected_results = { - 'AutoARIMA': 4.87, + 'AutoARIMA': 4.57, 'CES': 4.85, 'AutoETS': 4.35, 'DynamicOptimizedTheta': 4.54, - 'StatisticalEnsemble': 4.173 + 'StatisticalEnsemble': 4.23, } expected_results = pd.Series(expected_results) pd.testing.assert_series_equal( diff --git a/nbs/src/arima.ipynb b/nbs/src/arima.ipynb index 5dc05ab92..7ff394fc6 100644 --- a/nbs/src/arima.ipynb +++ b/nbs/src/arima.ipynb @@ -1,5 +1,17 @@ { "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "1f5499a6-5e03-48f4-aefd-226af8cffd62", + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, { "cell_type": "code", "execution_count": null, @@ -29,7 +41,7 @@ "outputs": [], "source": [ "#| hide\n", - "warnings.simplefilter('ignore')\n" + "warnings.simplefilter('ignore')" ] }, { @@ -57,12 +69,12 @@ "import numpy as np\n", "import pandas as pd\n", "import statsmodels.api as sm\n", - "from numba import njit\n", "from scipy.optimize import minimize\n", + "from scipy.signal import convolve\n", "from scipy.stats import norm\n", "\n", - "from statsforecast.mstl import mstl\n", - "from statsforecast.utils import CACHE, NOGIL" + "from statsforecast._lib import arima as _arima\n", + "from statsforecast.mstl import mstl" ] }, { @@ -91,75 +103,24 @@ { "cell_type": "code", "execution_count": null, - "id": "376c09fe", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def partrans(p, raw, new):\n", - " if p > 100:\n", - " raise ValueError('can only transform 100 pars in arima0')\n", - " \n", - " new[:p] = np.tanh(raw[:p])\n", - " work = new[:p].copy()\n", - " \n", - " for j in range(1, p):\n", - " a = new[j]\n", - " for k in range(j):\n", - " work[k] -= a * new[j - k - 1]\n", - " new[:j] = work[:j]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1ea975b7", + "id": "1d3b1233-779e-468f-8cd4-73f0ea54932e", "metadata": {}, "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def arima_gradtrans(x, arma):\n", - " eps = 1e-3\n", - " mp, mq, msp = arma[:3]\n", - " n = len(x)\n", - " y = np.identity(n)\n", - " w1 = np.empty(100)\n", - " w2 = np.empty(100)\n", - " w3 = np.empty(100)\n", - " if mp > 0:\n", - " for i in range(mp):\n", - " w1[i] = x[i]\n", - " partrans(mp, w1, w2)\n", - " for i in range(mp):\n", - " w1[i] += eps\n", - " partrans(mp, w1, w3)\n", - " for j in range(mp):\n", - " y[i, j] = (w3[j] - w2[j]) / eps\n", - " w1[i] -= eps\n", - " if msp > 0:\n", - " v = mp + mq\n", - " for i in range(msp):\n", - " w1[i] = x[i + v]\n", - " partrans(msp, w1, w2)\n", - " for j in range(msp):\n", - " w1[i] += eps\n", - " partrans(msp, w1, w3)\n", - " y[i + v, j + v] = (w3[j] - w2[j]) / eps\n", - " w1[i] -= eps\n", - " return y" + " return _arima.arima_gradtrans(x, arma)" ] }, { "cell_type": "code", "execution_count": null, - "id": "8c4d8301", + "id": "0d0b1c6d-5c05-409c-9d5e-81ea96cd8467", "metadata": {}, "outputs": [], "source": [ "#| hide\n", - "x = np.array([0.1, 0.4, 1.0, 3.1])\n", + "x = np.array([0.1, 0.4, 1.0, 3.1], dtype=np.float32)\n", "arma = np.array([1, 0, 1])\n", "expected = np.diag([0.9899673, 0.8553135, 1, 1])\n", "np.testing.assert_allclose(arima_gradtrans(x, arma), expected)" @@ -173,16 +134,8 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def arima_undopars(x, arma):\n", - " mp, mq, msp = arma[:3]\n", - " res = x.copy()\n", - " if mp > 0:\n", - " partrans(mp, x, res)\n", - " v = mp + mq\n", - " if msp > 0:\n", - " partrans(msp, x[v:], res[v:])\n", - " return res" + " return _arima.arima_undopars(x, arma)" ] }, { @@ -193,113 +146,12 @@ "outputs": [], "source": [ "#| hide\n", + "x = np.array([0.1, 0.4, 1.0, 3.1])\n", + "arma = np.array([1, 0, 1])\n", "expected = np.array([0.09966799, 0.37994896, 1.00000000, 3.10000000])\n", "np.testing.assert_allclose(arima_undopars(x, arma), expected)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "8ca693ea", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def tsconv(a, b):\n", - " na = len(a)\n", - " nb = len(b)\n", - " \n", - " nab = na + nb - 1\n", - " ab = np.zeros(nab)\n", - " \n", - " for i in range(na):\n", - " for j in range(nb):\n", - " ab[i + j] += a[i] * b[j]\n", - " \n", - " return ab" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a94fe482", - "metadata": {}, - "outputs": [], - "source": [ - "#| hide\n", - "x = np.arange(1, 11)\n", - "expected_tsconv = np.array([\n", - " 1, 4, 10, 20, 35, 56, 84, 120, 165, 220, 264,\n", - " 296, 315, 320, 310, 284, 241, 180, 100\n", - "])\n", - "\n", - "np.testing.assert_allclose(expected_tsconv, tsconv(x, x))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0f18a037", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def inclu2(np_, xnext, xrow, ynext, d, rbar, thetab):\n", - " for i in range(np_):\n", - " xrow[i] = xnext[i]\n", - " \n", - " ithisr = 0\n", - " for i in range(np_):\n", - " if xrow[i] != 0.:\n", - " xi = xrow[i]\n", - " di = d[i]\n", - " dpi = di + xi * xi\n", - " d[i] = dpi\n", - " cbar = di / dpi if dpi != 0. else math.inf\n", - " sbar = xi / dpi if dpi != 0. else math.inf\n", - " for k in range(i + 1, np_):\n", - " xk = xrow[k]\n", - " rbthis = rbar[ithisr]\n", - " xrow[k] = xk - xi * rbthis\n", - " rbar[ithisr] = cbar * rbthis + sbar * xk\n", - " ithisr += 1\n", - " xk = ynext\n", - " ynext = xk - xi * thetab[i]\n", - " thetab[i] = cbar * thetab[i] + sbar * xk\n", - " if di == 0.:\n", - " return\n", - " else:\n", - " ithisr = ithisr + np_ - i - 1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0db3588b", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def invpartrans(p, phi, new):\n", - " if p > 100:\n", - " raise ValueError('can only transform 100 pars in arima0')\n", - "\n", - " new = phi[:p].copy()\n", - " work = new.copy()\n", - " for k in range(p-1):\n", - " j = p - k - 1\n", - " a = new[j]\n", - " for k in range(j):\n", - " work[k] = (new[k] + a * new[j - k - 1]) / (1 - a * a)\n", - " for k in range(j):\n", - " new[k] = work[k]\n", - " for j in range(p):\n", - " new[j] = math.atanh(new[j])" - ] - }, { "cell_type": "code", "execution_count": null, @@ -308,15 +160,14 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def ARIMA_invtrans(x, arma):\n", " mp, mq, msp = arma[:3]\n", " y = x.copy()\n", " if mp > 0:\n", - " invpartrans(mp, x, y)\n", + " _arima.invpartrans(mp, x, y)\n", " v = mp + mq\n", " if msp > 0:\n", - " invpartrans(msp, x[v:], y[v:])\n", + " _arima.invpartrans(msp, x[v:], y[v:])\n", " return y" ] }, @@ -341,134 +192,13 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def getQ0(phi, theta):\n", " p = len(phi)\n", " q = len(theta)\n", " r = max(p, q + 1)\n", - " \n", - " np_ = r * (r + 1) // 2\n", - " nrbar = np_ * (np_ - 1) // 2\n", - " \n", - " V = np.zeros(np_)\n", - " ind = 0\n", - " for j in range(r):\n", - " vj = 0.\n", - " if j == 0:\n", - " vj = 1.\n", - " elif j - 1 < q:\n", - " vj = theta[j - 1]\n", - " \n", - " for i in range(j, r):\n", - " vi = 0.\n", - " if i == 0:\n", - " vi = 1.0\n", - " elif i - 1 < q:\n", - " vi = theta[i - 1]\n", - " V[ind] = vi * vj\n", - " ind += 1\n", - " \n", - " res = np.zeros((r, r))\n", - " res = res.flatten()\n", - " \n", - " if r == 1:\n", - " if p == 0:\n", - " res[0] = 1.\n", - " else:\n", - " res[0] = 1. / (1. - phi[0] * phi[0])\n", - " \n", - " res = res.reshape((r, r))\n", - " return res\n", - " \n", - " if p > 0:\n", - " rbar = np.zeros(nrbar)\n", - " thetab = np.zeros(np_)\n", - " xnext = np.zeros(np_)\n", - " xrow = np.zeros(np_)\n", - " \n", - " ind = 0\n", - " ind1 = -1\n", - " npr = np_ - r\n", - " npr1 = npr + 1\n", - " indj = npr\n", - " ind2 = npr - 1\n", - " \n", - " for j in range(r):\n", - " phij = phi[j] if j < p else 0.\n", - " xnext[indj] = 0.\n", - " indj += 1\n", - " indi = npr1 + j\n", - " for i in range(j, r):\n", - " ynext = V[ind]\n", - " ind += 1\n", - " phii = phi[i] if i < p else 0.\n", - " if j != r - 1:\n", - " xnext[indj] = -phii\n", - " if i != r - 1:\n", - " xnext[indi] -= phij\n", - " ind1 += 1\n", - " xnext[ind1] = -1.\n", - " xnext[npr] = -phii * phij\n", - " ind2 += 1\n", - " if ind2 >= np_:\n", - " ind2 = 0\n", - " xnext[ind2] += 1.\n", - " inclu2(np_, xnext, xrow, ynext, res, rbar, thetab)\n", - " xnext[ind2] = 0.\n", - " if i != r - 1:\n", - " xnext[indi] = 0.\n", - " indi += 1\n", - " xnext[ind1] = 0.\n", - " \n", - " ithisr = nrbar - 1\n", - " im = np_ - 1\n", - " for i in range(np_):\n", - " bi = thetab[im]\n", - " jm = np_ - 1\n", - " for j in range(i):\n", - " bi -= rbar[ithisr] * res[jm]\n", - " ithisr -= 1\n", - " jm -= 1\n", - " res[im] = bi\n", - " im -= 1\n", - " \n", - " # Now reorder p\n", - " ind = npr\n", - " for i in range(r):\n", - " xnext[i] = res[ind]\n", - " ind += 1\n", - " ind = np_ - 1\n", - " ind1 = npr - 1\n", - " for i in range(npr):\n", - " res[ind] = res[ind1]\n", - " ind -= 1\n", - " ind1 -= 1\n", - " for i in range(r):\n", - " res[i] = xnext[i]\n", - " else:\n", - " indn = np_\n", - " ind = np_\n", - " for i in range(r):\n", - " for j in range(i + 1):\n", - " ind -= 1\n", - " res[ind] = V[ind]\n", - " if j != 0:\n", - " indn -= 1\n", - " res[ind] += res[ind]\n", - " \n", - " # Unpack to a full matrix\n", - " ind = np_\n", - " for i in range(r - 1, 0, -1):\n", - " for j in range(r - 1, i - 1, -1):\n", - " ind -= 1\n", - " res[r * i + j] = res[ind]\n", - "\n", - " for i in range(r - 1):\n", - " for j in range(i + 1, r):\n", - " res[i + r * j] = res[j + r * i]\n", - " \n", - " res = res.reshape((r, r))\n", - " return res" + " res = np.zeros(r * r, dtype=np.float64)\n", + " _arima.getQ0(phi, theta, res)\n", + " return res.reshape(r, r)" ] }, { @@ -526,43 +256,9 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def arima_transpar(params_in, arma, trans):\n", " #TODO check trans=True results\n", - " mp, mq, msp, msq, ns = arma[:5]\n", - " p = mp + ns * msp\n", - " q = mq + ns * msq\n", - " \n", - " phi = np.zeros(p)\n", - " theta = np.zeros(q)\n", - " params = params_in.copy()\n", - " \n", - " if trans:\n", - " #n = mp + mq + msp + msq\n", - " if mp > 0:\n", - " partrans(mp, params_in, params)\n", - " v = mp + mq\n", - " if msp > 0:\n", - " partrans(msp, params_in[v:], params[v:])\n", - " if ns > 0:\n", - " phi[:mp] = params[:mp]\n", - " phi[mp:p] = 0.\n", - " theta[:mq] = params[mp:mp+mq]\n", - " theta[mq:q] = 0.\n", - " for j in range(msp):\n", - " phi[(j + 1) * ns - 1] += params[j + mp + mq]\n", - " for i in range(mp):\n", - " phi[(j + 1) * ns + i] -= params[i] * params[j + mp + mq]\n", - " \n", - " for j in range(msq):\n", - " theta[(j + 1) * ns - 1] += params[j + mp + mq + msp]\n", - " for i in range(mq):\n", - " theta[(j + 1) * ns + i] += params[i + mp] * params[j + mp + mq + msp]\n", - " else:\n", - " phi[:mp] = params[:mp]\n", - " theta[:mq] = theta[mp:mp + mq]\n", - " \n", - " return phi, theta" + " return _arima.arima_transpar(params_in, arma, trans)" ] }, { @@ -608,48 +304,8 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def arima_css(y, arma, phi, theta, ncond):\n", - " n = len(y)\n", - " p = len(phi)\n", - " q = len(theta)\n", - " nu = 0\n", - " ssq = 0.0\n", - " \n", - " w = y.copy()\n", - " \n", - " for i in range(arma[5]):\n", - " for l in range(n - 1, 0, -1):\n", - " w[l] -= w[l - 1]\n", - " \n", - " ns = arma[4]\n", - " for i in range(arma[6]):\n", - " for l in range(n - 1, ns - 1, -1):\n", - " w[l] -= w[l - ns]\n", - " \n", - " resid = np.empty(n)\n", - " resid[:ncond] = 0.\n", - " for l in range(ncond, n):\n", - " tmp = w[l]\n", - " for j in range(p):\n", - " if l - j - 1 < 0:\n", - " continue\n", - " tmp -= phi[j] * w[l - j - 1]\n", - " \n", - " for j in range(min(l - ncond, q)):\n", - " if l - j - 1 < 0:\n", - " continue\n", - " tmp -= theta[j] * resid[l - j - 1]\n", - " \n", - " resid[l] = tmp\n", - " \n", - " if not np.isnan(tmp):\n", - " nu += 1\n", - " ssq += tmp * tmp\n", - " \n", - " res = ssq / nu\n", - " \n", - " return res, resid" + "def arima_css(y, arma, phi, theta):\n", + " return _arima.arima_css(y, arma, phi, theta)" ] }, { @@ -660,11 +316,12 @@ "outputs": [], "source": [ "#| hide\n", - "arima_css(np.arange(1, 11), \n", - " np.array([0,0,0,0,0,0,0], dtype=np.int32),\n", - " expected_arima_transpar_f[0],\n", - " expected_arima_transpar_f[1], \n", - " 3)" + "arima_css(\n", + " np.arange(1, 11), \n", + " np.array([0,0,0,0,0,0,0], dtype=np.int32),\n", + " expected_arima_transpar_f[0],\n", + " expected_arima_transpar_f[1]\n", + ")" ] }, { @@ -675,8 +332,7 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def _make_arima(phi, theta, delta, kappa = 1e6, tol = np.finfo(float).eps):\n", + "def make_arima(phi, theta, delta, kappa = 1e6, tol = np.finfo(float).eps):\n", " # check nas phi\n", " # check nas theta\n", " p = len(phi)\n", @@ -691,14 +347,14 @@ " if p > 0:\n", " T[:p, 0] = phi\n", " if r > 1:\n", - " for i in range(1, r):\n", - " T[i - 1, i] = 1\n", + " idx = np.arange(1, r)\n", + " T[idx - 1, idx] = 1\n", "\n", " if d > 0:\n", " T[r] = Z\n", " if d > 1:\n", - " for ind in range(1, d):\n", - " T[r + ind, r + ind - 1] = 1\n", + " idx = np.arange(1, d)\n", + " T[r + idx, r + idx - 1] = 1\n", "\n", " if q < r - 1:\n", " theta = np.concatenate((theta, np.zeros(r - 1 - q)))\n", @@ -716,15 +372,21 @@ " Pn[0, 0] = 1 / (1 - phi[0] ** 2) if p > 0 else 1.\n", " \n", " if d > 0:\n", - " for i in range(d):\n", - " Pn[r + i, r + i] = kappa\n", + " idx = np.arange(d)\n", + " Pn[r + idx, r + idx] = kappa\n", " \n", - " return phi, theta, delta, Z, a, P, T, V, h, Pn\n", - "\n", - "def make_arima(phi, theta, delta, kappa = 1e6, tol = np.finfo(np.float64).eps):\n", - " keys = ['phi', 'theta', 'delta', 'Z', 'a', 'P', 'T', 'V', 'h', 'Pn']\n", - " res = _make_arima(phi, theta, delta, kappa, tol)\n", - " return dict(zip(keys, res))" + " return {\n", + " 'phi': phi,\n", + " 'theta': theta,\n", + " 'delta': delta,\n", + " 'Z': Z,\n", + " 'a': a,\n", + " 'P': P,\n", + " 'T': T,\n", + " 'V': V,\n", + " 'h': h,\n", + " 'Pn': Pn,\n", + " }" ] }, { @@ -749,134 +411,23 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def arima_like(y, phi, theta, delta, a, P, Pn, up, use_resid):\n", - " n = len(y)\n", - " rd = len(a)\n", - " p = len(phi)\n", - " q = len(theta)\n", - " d = len(delta)\n", - " r = rd - d\n", - " \n", - " sumlog = 0.\n", - " ssq = 0.\n", - " nu = 0\n", - " \n", - " P = P.ravel()\n", - " Pnew = Pn.ravel()\n", - " anew = np.empty(rd)\n", - " M = np.empty(rd)\n", - " if d > 0:\n", - " mm = np.empty(rd * rd)\n", - "\n", " if use_resid:\n", - " rsResid = np.empty(n)\n", - " \n", - " for l in range(n):\n", - " for i in range(r):\n", - " tmp = a[i + 1] if i < r - 1 else 0.\n", - " if i < p:\n", - " tmp += phi[i] * a[0]\n", - " anew[i] = tmp\n", - " if d > 0:\n", - " for i in range(r + 1, rd):\n", - " anew[i] = a[i - 1]\n", - " tmp = a[0]\n", - " for i in range(d):\n", - " tmp += delta[i] * a[r + i]\n", - " anew[r] = tmp\n", - " if l > up:\n", - " if d == 0:\n", - " for i in range(r):\n", - " vi = 0.\n", - " if i == 0:\n", - " vi = 1.\n", - " elif i - 1 < q:\n", - " vi = theta[i - 1]\n", - " for j in range(r):\n", - " tmp = 0.\n", - " if j == 0:\n", - " tmp = vi\n", - " elif j - 1 < q:\n", - " tmp = vi * theta[j - 1]\n", - " if i < p and j < p:\n", - " tmp += phi[i] * phi[j] * P[0]\n", - " if i < r - 1 and j < r -1:\n", - " tmp += P[i + 1 + r * (j + 1)]\n", - " if i < p and j < r - 1:\n", - " tmp += phi[i] * P[j + 1]\n", - " if j < p and i < r -1:\n", - " tmp += phi[j] * P[i + 1]\n", - " Pnew[i + r * j] = tmp\n", - " else:\n", - " # mm = TP\n", - " for i in range(r):\n", - " for j in range(rd):\n", - " tmp = 0.\n", - " if i < p:\n", - " tmp += phi[i] * P[rd * j]\n", - " if i < r - 1:\n", - " tmp += P[i + 1 + rd * j]\n", - " mm[i + rd * j] = tmp\n", - " for j in range(rd):\n", - " tmp = P[rd * j]\n", - " for k in range(d):\n", - " tmp += delta[k] * P[r + k + rd * j]\n", - " mm[r + rd * j] = tmp\n", - " for i in range(1, d):\n", - " for j in range(rd):\n", - " mm[r + i + rd * j] = P[r + i - 1 + rd * j]\n", - " \n", - " # Pnew = mmT'\n", - " for i in range(r):\n", - " for j in range(rd):\n", - " tmp = 0.\n", - " if i < p:\n", - " tmp += phi[i] * mm[j]\n", - " if i < r - 1:\n", - " tmp += mm[rd * (i + 1) + j]\n", - " Pnew[j + rd * i] = tmp\n", - " for j in range(rd):\n", - " tmp = mm[j]\n", - " for k in range(d):\n", - " tmp += delta[k] * mm[rd * (r + k) + j]\n", - " Pnew[rd * r + j] = tmp\n", - " for i in range(1, d):\n", - " for j in range(rd):\n", - " Pnew[rd * (r + i) + j] = mm[rd * (r + i - 1) + j]\n", - " for i in range(q + 1):\n", - " vi = 1. if i == 0 else theta[i - 1]\n", - " for j in range(q + 1):\n", - " Pnew[i + rd * j] += vi * (1. if j == 0 else theta[j - 1])\n", - " \n", - " if not math.isnan(y[l]):\n", - " resid = y[l] - anew[0]\n", - " for i in range(d):\n", - " resid -= delta[i] * anew[r + i]\n", - " for i in range(rd):\n", - " tmp = Pnew[i]\n", - " for j in range(d):\n", - " tmp += Pnew[i + (r + j) * rd] * delta[j]\n", - " M[i] = tmp\n", - " gain = M[0]\n", - " for j in range(d):\n", - " gain += delta[j] * M[r + j]\n", - " if gain < 1e4:\n", - " nu += 1\n", - " ssq += resid * resid / gain if gain != 0. else math.inf\n", - " sumlog += math.log(gain)\n", - " if use_resid:\n", - " rsResid[l] = resid / math.sqrt(gain) if gain != 0. else math.inf\n", - " for i in range(rd):\n", - " a[i] = anew[i] + M[i] * resid / gain if gain != 0. else math.inf\n", - " for i in range(rd):\n", - " for j in range(rd):\n", - " P[i + j * rd] = Pnew[i + j * rd] - M[i] * M[j] / gain if gain != 0. else math.inf\n", - " else:\n", - " a[:] = anew[:]\n", - " P[:] = Pnew[:]\n", - " if use_resid:\n", - " rsResid[l] = np.nan\n", + " rsResid = np.empty_like(y)\n", + " else:\n", + " rsResid = np.empty(0)\n", + " ssq, sumlog, nu = _arima.arima_like(\n", + " y,\n", + " phi,\n", + " theta,\n", + " delta,\n", + " a,\n", + " P.ravel(),\n", + " Pn.ravel(),\n", + " up,\n", + " use_resid,\n", + " rsResid,\n", + " )\n", " if not use_resid:\n", " rsResid = None\n", " return ssq, sumlog, nu, rsResid" @@ -890,7 +441,7 @@ "outputs": [], "source": [ "#| hide\n", - "y = np.arange(10)\n", + "y = np.arange(10, dtype=np.float64)\n", "phi = np.array([0.99551517])\n", "theta = np.array([])\n", "delta = np.array([1.0])\n", @@ -902,7 +453,8 @@ "])\n", "up = 0\n", "use_resid = True\n", - "res = arima_like(y, phi, theta, delta, a, P, Pn, up, use_resid)" + "res = arima_like(y, phi, theta, delta, a, P, Pn, up, use_resid)\n", + "res" ] }, { @@ -913,35 +465,15 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def diff1d(x, lag, differences):\n", + "def diff(x, lag, differences):\n", + " x = np.asarray(x, dtype=np.float64)\n", " y = x.copy()\n", " for _ in range(differences):\n", " x = y.copy()\n", - " for i in range(lag):\n", - " y[i] = np.nan\n", - " for i in range(lag, x.size):\n", - " y[i] = x[i] - x[i - lag]\n", - " return y\n", - "\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def diff2d(x, lag, differences):\n", - " y = np.empty_like(x)\n", - " for j in range(x.shape[1]):\n", - " y[:, j] = diff1d(x[:, j], lag, differences)\n", - " return y\n", - "\n", - "\n", - "def diff(x, lag, differences):\n", - " if x.ndim == 1:\n", - " y = diff1d(x, lag, differences)\n", - " nan_mask = np.isnan(y)\n", - " elif x.ndim == 2:\n", - " y = diff2d(x, lag, differences)\n", - " nan_mask = np.isnan(y).all(1)\n", - " else:\n", - " raise ValueError(x.ndim)\n", - " return y[~nan_mask]" + " y[:lag] = np.nan\n", + " y[lag:] = x[lag:] - x[:-lag]\n", + " nans = lag * differences\n", + " return y[nans:]" ] }, { @@ -1009,7 +541,7 @@ " tol=1e-8,\n", " optim_control = {'maxiter': 100}):\n", " SSG = SSinit == 'Gardner1980'\n", - " x = x.copy()\n", + " x = x.astype(np.float64, copy=True)\n", " \n", " def upARIMA(mod, phi, theta):\n", " p = len(phi)\n", @@ -1118,21 +650,26 @@ " \n", " #fixed\n", " #mask \n", - " arma = (*order[::2], \n", + " arma = np.array(\n", + " [\n", + " *order[::2], \n", " *seasonal['order'][::2],\n", " seasonal['period'],\n", " order[1],\n", - " seasonal['order'][1])\n", - " narma = sum(arma[:4])\n", + " seasonal['order'][1],\n", + " ],\n", + " dtype=np.intc,\n", + " )\n", + " narma = arma[:4].sum().item()\n", " \n", " # xtsp = init x, end x and frequency\n", " # tsp(x) = None\n", " Delta = np.array([1.]) \n", " for i in range(order[1]):\n", - " Delta = tsconv(Delta, np.array([1., -1.])) \n", + " Delta = convolve(Delta, np.array([1., -1.])) \n", " \n", " for i in range(seasonal['order'][1]):\n", - " Delta = tsconv(Delta, np.array([1] + [0]*(seasonal['period'] - 1) + [-1]))\n", + " Delta = convolve(Delta, np.array([1] + [0]*(seasonal['period'] - 1) + [-1]))\n", " Delta = - Delta[1:]\n", " nd = order[1] + seasonal['order'][1]\n", " n_used = (~np.isnan(x)).sum() - len(Delta)\n", @@ -1257,7 +794,7 @@ " if ncxreg > 0:\n", " x -= np.dot(xreg, par[narma + np.arange(ncxreg)])\n", "\n", - " res, resid = arima_css(x, arma, phi, theta, ncond)\n", + " res, _ = arima_css(x, arma, phi, theta)\n", " if math.isinf(res):\n", " import sys\n", "\n", @@ -1284,7 +821,7 @@ " mod = make_arima(phi, theta, Delta, kappa)\n", " if ncxreg > 0:\n", " x -= np.dot(xreg, coef[narma + np.arange(ncxreg)])\n", - " val = arima_css(x, arma, phi, theta, ncond)\n", + " val = arima_css(x, arma, phi, theta)\n", " sigma2 = val[0]\n", " var = None if no_optim else res.hess_inv / n_used\n", " else:\n", @@ -1295,7 +832,7 @@ " # only update the initial parameters if they're valid\n", " candidate = init.copy()\n", " candidate[mask] = res.x\n", - " phi, _ = arima_transpar(candidate, arma, False)\n", + " phi, _ = arima_transpar(candidate, arma, transform_pars)\n", " if np.logical_and(phi > - math.pi / 2, phi < math.pi / 2).all():\n", " init = candidate\n", " if arma[0] > 0:\n", @@ -1306,13 +843,16 @@ " raise ValueError('non-stationary seasonal AR part from CSS')\n", " ncond = 0\n", " if transform_pars:\n", - " init = ARIMA_invtrans(init, arma)\n", - " if arma[1] > 0:\n", - " ind = arma[0] + np.arange(arma[1])\n", - " init[ind] = maInvert(init[ind])\n", - " if arma[3] > 0:\n", - " ind = np.sum(arma[:3]) + np.arange(arma[3])\n", - " init[ind] = maInvert(init[ind])\n", + " candidate = ARIMA_invtrans(init, arma)\n", + " phi, _ = arima_transpar(candidate, arma, transform_pars)\n", + " if np.logical_and(phi > - math.pi / 2, phi < math.pi / 2).all():\n", + " init = candidate\n", + " if arma[1] > 0:\n", + " ind = arma[0] + np.arange(arma[1])\n", + " init[ind] = maInvert(init[ind])\n", + " if arma[3] > 0:\n", + " ind = np.sum(arma[:3]) + np.arange(arma[3])\n", + " init[ind] = maInvert(init[ind])\n", " trarma = arima_transpar(init, arma, transform_pars)\n", " mod = make_arima(trarma[0], trarma[1], Delta, kappa, SSinit)\n", " if no_optim:\n", @@ -1379,7 +919,7 @@ " 'mask': mask,\n", " 'loglik': -0.5 * value, \n", " 'aic': aic, \n", - " 'arma': arma,\n", + " 'arma': tuple(x.item() for x in arma),\n", " 'residuals': resid, \n", " #'series': series,\n", " 'code': res.status, \n", @@ -1511,46 +1051,19 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def kalman_forecast(n, Z, a, P, T, V, h):\n", - " p = len(a)\n", - " \n", " a = a.copy()\n", - " anew = np.empty(p)\n", - " Pnew = np.empty((p, p))\n", - " mm = np.empty((p, p))\n", + " P = P.copy() \n", " forecasts = np.empty(n)\n", " se = np.empty(n)\n", - " P = P.copy()\n", - " \n", + " z = Z.reshape(-1, 1) * Z.reshape(1, -1)\n", " for l in range(n):\n", - " anew = T @ a\n", - " \n", - " a[:] = anew[:]\n", - " forecasts[l] = anew @ Z\n", - " \n", - " for i in range(p):\n", - " for j in range(p):\n", - " tmp = 0.\n", - " for k in range(p):\n", - " tmp += T[i, k] * P[k, j]\n", - " mm[i, j] = tmp\n", - "\n", - " for i in range(p):\n", - " for j in range(p):\n", - " tmp = V[i, j]\n", - " for k in range(p):\n", - " tmp += mm[i, k] * T[j, k]\n", - " Pnew[i, j] = tmp\n", - "\n", - " tmp = h\n", - " for i in range(p):\n", - " for j in range(p):\n", - " P[i, j] = Pnew[i, j]\n", - " tmp += Z[i] * Z[j] * P[i, j]\n", - " se[l] = tmp\n", - "\n", - " return forecasts, se" + " a = T @ a\n", + " forecasts[l] = a @ Z\n", + " mm = T @ P\n", + " P = V + mm @ T.T\n", + " se[l] = h + np.sum(z * P)\n", + " return forecasts, se " ] }, { diff --git a/python/statsforecast/_modidx.py b/python/statsforecast/_modidx.py index 246699083..2a6176659 100644 --- a/python/statsforecast/_modidx.py +++ b/python/statsforecast/_modidx.py @@ -32,7 +32,6 @@ 'statsforecast/arima.py'), 'statsforecast.arima.AutoARIMA.summary': ( 'src/arima.html#autoarima.summary', 'statsforecast/arima.py'), - 'statsforecast.arima._make_arima': ('src/arima.html#_make_arima', 'statsforecast/arima.py'), 'statsforecast.arima.arima': ('src/arima.html#arima', 'statsforecast/arima.py'), 'statsforecast.arima.arima2': ('src/arima.html#arima2', 'statsforecast/arima.py'), 'statsforecast.arima.arima_css': ('src/arima.html#arima_css', 'statsforecast/arima.py'), @@ -48,16 +47,12 @@ 'statsforecast.arima.convert_coef_name': ( 'src/arima.html#convert_coef_name', 'statsforecast/arima.py'), 'statsforecast.arima.diff': ('src/arima.html#diff', 'statsforecast/arima.py'), - 'statsforecast.arima.diff1d': ('src/arima.html#diff1d', 'statsforecast/arima.py'), - 'statsforecast.arima.diff2d': ('src/arima.html#diff2d', 'statsforecast/arima.py'), 'statsforecast.arima.fitted_arima': ('src/arima.html#fitted_arima', 'statsforecast/arima.py'), 'statsforecast.arima.fixed_params_from_dict': ( 'src/arima.html#fixed_params_from_dict', 'statsforecast/arima.py'), 'statsforecast.arima.forecast_arima': ('src/arima.html#forecast_arima', 'statsforecast/arima.py'), 'statsforecast.arima.forward_arima': ('src/arima.html#forward_arima', 'statsforecast/arima.py'), 'statsforecast.arima.getQ0': ('src/arima.html#getq0', 'statsforecast/arima.py'), - 'statsforecast.arima.inclu2': ('src/arima.html#inclu2', 'statsforecast/arima.py'), - 'statsforecast.arima.invpartrans': ('src/arima.html#invpartrans', 'statsforecast/arima.py'), 'statsforecast.arima.is_constant': ('src/arima.html#is_constant', 'statsforecast/arima.py'), 'statsforecast.arima.kalman_forecast': ('src/arima.html#kalman_forecast', 'statsforecast/arima.py'), 'statsforecast.arima.make_arima': ('src/arima.html#make_arima', 'statsforecast/arima.py'), @@ -65,13 +60,11 @@ 'statsforecast.arima.ndiffs': ('src/arima.html#ndiffs', 'statsforecast/arima.py'), 'statsforecast.arima.newmodel': ('src/arima.html#newmodel', 'statsforecast/arima.py'), 'statsforecast.arima.nsdiffs': ('src/arima.html#nsdiffs', 'statsforecast/arima.py'), - 'statsforecast.arima.partrans': ('src/arima.html#partrans', 'statsforecast/arima.py'), 'statsforecast.arima.predict_arima': ('src/arima.html#predict_arima', 'statsforecast/arima.py'), 'statsforecast.arima.print_statsforecast_ARIMA': ( 'src/arima.html#print_statsforecast_arima', 'statsforecast/arima.py'), 'statsforecast.arima.search_arima': ('src/arima.html#search_arima', 'statsforecast/arima.py'), - 'statsforecast.arima.seas_heuristic': ('src/arima.html#seas_heuristic', 'statsforecast/arima.py'), - 'statsforecast.arima.tsconv': ('src/arima.html#tsconv', 'statsforecast/arima.py')}, + 'statsforecast.arima.seas_heuristic': ('src/arima.html#seas_heuristic', 'statsforecast/arima.py')}, 'statsforecast.ces': { 'statsforecast.ces._simulate_pred_intervals': ( 'src/ces.html#_simulate_pred_intervals', 'statsforecast/ces.py'), 'statsforecast.ces.auto_ces': ('src/ces.html#auto_ces', 'statsforecast/ces.py'), diff --git a/python/statsforecast/arima.py b/python/statsforecast/arima.py index 42c9b5255..8cc7ebb61 100644 --- a/python/statsforecast/arima.py +++ b/python/statsforecast/arima.py @@ -4,7 +4,7 @@ __all__ = ['predict_arima', 'arima_string', 'forecast_arima', 'fitted_arima', 'auto_arima_f', 'print_statsforecast_ARIMA', 'ARIMASummary', 'AutoARIMA'] -# %% ../../nbs/src/arima.ipynb 4 +# %% ../../nbs/src/arima.ipynb 5 import math import warnings from collections import namedtuple @@ -14,365 +14,55 @@ import numpy as np import pandas as pd import statsmodels.api as sm -from numba import njit from scipy.optimize import minimize +from scipy.signal import convolve from scipy.stats import norm +from ._lib import arima as _arima from .mstl import mstl -from .utils import CACHE, NOGIL - -# %% ../../nbs/src/arima.ipynb 6 -OptimResult = namedtuple("OptimResult", "success status x fun hess_inv") # %% ../../nbs/src/arima.ipynb 7 -@njit(nogil=NOGIL, cache=CACHE) -def partrans(p, raw, new): - if p > 100: - raise ValueError("can only transform 100 pars in arima0") - - new[:p] = np.tanh(raw[:p]) - work = new[:p].copy() - - for j in range(1, p): - a = new[j] - for k in range(j): - work[k] -= a * new[j - k - 1] - new[:j] = work[:j] +OptimResult = namedtuple("OptimResult", "success status x fun hess_inv") # %% ../../nbs/src/arima.ipynb 8 -@njit(nogil=NOGIL, cache=CACHE) def arima_gradtrans(x, arma): - eps = 1e-3 - mp, mq, msp = arma[:3] - n = len(x) - y = np.identity(n) - w1 = np.empty(100) - w2 = np.empty(100) - w3 = np.empty(100) - if mp > 0: - for i in range(mp): - w1[i] = x[i] - partrans(mp, w1, w2) - for i in range(mp): - w1[i] += eps - partrans(mp, w1, w3) - for j in range(mp): - y[i, j] = (w3[j] - w2[j]) / eps - w1[i] -= eps - if msp > 0: - v = mp + mq - for i in range(msp): - w1[i] = x[i + v] - partrans(msp, w1, w2) - for j in range(msp): - w1[i] += eps - partrans(msp, w1, w3) - y[i + v, j + v] = (w3[j] - w2[j]) / eps - w1[i] -= eps - return y + return _arima.arima_gradtrans(x, arma) # %% ../../nbs/src/arima.ipynb 10 -@njit(nogil=NOGIL, cache=CACHE) def arima_undopars(x, arma): - mp, mq, msp = arma[:3] - res = x.copy() - if mp > 0: - partrans(mp, x, res) - v = mp + mq - if msp > 0: - partrans(msp, x[v:], res[v:]) - return res + return _arima.arima_undopars(x, arma) # %% ../../nbs/src/arima.ipynb 12 -@njit(nogil=NOGIL, cache=CACHE) -def tsconv(a, b): - na = len(a) - nb = len(b) - - nab = na + nb - 1 - ab = np.zeros(nab) - - for i in range(na): - for j in range(nb): - ab[i + j] += a[i] * b[j] - - return ab - -# %% ../../nbs/src/arima.ipynb 14 -@njit(nogil=NOGIL, cache=CACHE) -def inclu2(np_, xnext, xrow, ynext, d, rbar, thetab): - for i in range(np_): - xrow[i] = xnext[i] - - ithisr = 0 - for i in range(np_): - if xrow[i] != 0.0: - xi = xrow[i] - di = d[i] - dpi = di + xi * xi - d[i] = dpi - cbar = di / dpi if dpi != 0.0 else math.inf - sbar = xi / dpi if dpi != 0.0 else math.inf - for k in range(i + 1, np_): - xk = xrow[k] - rbthis = rbar[ithisr] - xrow[k] = xk - xi * rbthis - rbar[ithisr] = cbar * rbthis + sbar * xk - ithisr += 1 - xk = ynext - ynext = xk - xi * thetab[i] - thetab[i] = cbar * thetab[i] + sbar * xk - if di == 0.0: - return - else: - ithisr = ithisr + np_ - i - 1 - -# %% ../../nbs/src/arima.ipynb 15 -@njit(nogil=NOGIL, cache=CACHE) -def invpartrans(p, phi, new): - if p > 100: - raise ValueError("can only transform 100 pars in arima0") - - new = phi[:p].copy() - work = new.copy() - for k in range(p - 1): - j = p - k - 1 - a = new[j] - for k in range(j): - work[k] = (new[k] + a * new[j - k - 1]) / (1 - a * a) - for k in range(j): - new[k] = work[k] - for j in range(p): - new[j] = math.atanh(new[j]) - -# %% ../../nbs/src/arima.ipynb 16 -@njit(nogil=NOGIL, cache=CACHE) def ARIMA_invtrans(x, arma): mp, mq, msp = arma[:3] y = x.copy() if mp > 0: - invpartrans(mp, x, y) + _arima.invpartrans(mp, x, y) v = mp + mq if msp > 0: - invpartrans(msp, x[v:], y[v:]) + _arima.invpartrans(msp, x[v:], y[v:]) return y -# %% ../../nbs/src/arima.ipynb 18 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/arima.ipynb 14 def getQ0(phi, theta): p = len(phi) q = len(theta) r = max(p, q + 1) + res = np.zeros(r * r, dtype=np.float64) + _arima.getQ0(phi, theta, res) + return res.reshape(r, r) - np_ = r * (r + 1) // 2 - nrbar = np_ * (np_ - 1) // 2 - - V = np.zeros(np_) - ind = 0 - for j in range(r): - vj = 0.0 - if j == 0: - vj = 1.0 - elif j - 1 < q: - vj = theta[j - 1] - - for i in range(j, r): - vi = 0.0 - if i == 0: - vi = 1.0 - elif i - 1 < q: - vi = theta[i - 1] - V[ind] = vi * vj - ind += 1 - - res = np.zeros((r, r)) - res = res.flatten() - - if r == 1: - if p == 0: - res[0] = 1.0 - else: - res[0] = 1.0 / (1.0 - phi[0] * phi[0]) - - res = res.reshape((r, r)) - return res - - if p > 0: - rbar = np.zeros(nrbar) - thetab = np.zeros(np_) - xnext = np.zeros(np_) - xrow = np.zeros(np_) - - ind = 0 - ind1 = -1 - npr = np_ - r - npr1 = npr + 1 - indj = npr - ind2 = npr - 1 - - for j in range(r): - phij = phi[j] if j < p else 0.0 - xnext[indj] = 0.0 - indj += 1 - indi = npr1 + j - for i in range(j, r): - ynext = V[ind] - ind += 1 - phii = phi[i] if i < p else 0.0 - if j != r - 1: - xnext[indj] = -phii - if i != r - 1: - xnext[indi] -= phij - ind1 += 1 - xnext[ind1] = -1.0 - xnext[npr] = -phii * phij - ind2 += 1 - if ind2 >= np_: - ind2 = 0 - xnext[ind2] += 1.0 - inclu2(np_, xnext, xrow, ynext, res, rbar, thetab) - xnext[ind2] = 0.0 - if i != r - 1: - xnext[indi] = 0.0 - indi += 1 - xnext[ind1] = 0.0 - - ithisr = nrbar - 1 - im = np_ - 1 - for i in range(np_): - bi = thetab[im] - jm = np_ - 1 - for j in range(i): - bi -= rbar[ithisr] * res[jm] - ithisr -= 1 - jm -= 1 - res[im] = bi - im -= 1 - - # Now reorder p - ind = npr - for i in range(r): - xnext[i] = res[ind] - ind += 1 - ind = np_ - 1 - ind1 = npr - 1 - for i in range(npr): - res[ind] = res[ind1] - ind -= 1 - ind1 -= 1 - for i in range(r): - res[i] = xnext[i] - else: - indn = np_ - ind = np_ - for i in range(r): - for j in range(i + 1): - ind -= 1 - res[ind] = V[ind] - if j != 0: - indn -= 1 - res[ind] += res[ind] - - # Unpack to a full matrix - ind = np_ - for i in range(r - 1, 0, -1): - for j in range(r - 1, i - 1, -1): - ind -= 1 - res[r * i + j] = res[ind] - - for i in range(r - 1): - for j in range(i + 1, r): - res[i + r * j] = res[j + r * i] - - res = res.reshape((r, r)) - return res - -# %% ../../nbs/src/arima.ipynb 20 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/arima.ipynb 16 def arima_transpar(params_in, arma, trans): # TODO check trans=True results - mp, mq, msp, msq, ns = arma[:5] - p = mp + ns * msp - q = mq + ns * msq - - phi = np.zeros(p) - theta = np.zeros(q) - params = params_in.copy() - - if trans: - # n = mp + mq + msp + msq - if mp > 0: - partrans(mp, params_in, params) - v = mp + mq - if msp > 0: - partrans(msp, params_in[v:], params[v:]) - if ns > 0: - phi[:mp] = params[:mp] - phi[mp:p] = 0.0 - theta[:mq] = params[mp : mp + mq] - theta[mq:q] = 0.0 - for j in range(msp): - phi[(j + 1) * ns - 1] += params[j + mp + mq] - for i in range(mp): - phi[(j + 1) * ns + i] -= params[i] * params[j + mp + mq] - - for j in range(msq): - theta[(j + 1) * ns - 1] += params[j + mp + mq + msp] - for i in range(mq): - theta[(j + 1) * ns + i] += params[i + mp] * params[j + mp + mq + msp] - else: - phi[:mp] = params[:mp] - theta[:mq] = theta[mp : mp + mq] - - return phi, theta - -# %% ../../nbs/src/arima.ipynb 23 -@njit(nogil=NOGIL, cache=CACHE) -def arima_css(y, arma, phi, theta, ncond): - n = len(y) - p = len(phi) - q = len(theta) - nu = 0 - ssq = 0.0 - - w = y.copy() - - for i in range(arma[5]): - for l in range(n - 1, 0, -1): - w[l] -= w[l - 1] - - ns = arma[4] - for i in range(arma[6]): - for l in range(n - 1, ns - 1, -1): - w[l] -= w[l - ns] - - resid = np.empty(n) - resid[:ncond] = 0.0 - for l in range(ncond, n): - tmp = w[l] - for j in range(p): - if l - j - 1 < 0: - continue - tmp -= phi[j] * w[l - j - 1] - - for j in range(min(l - ncond, q)): - if l - j - 1 < 0: - continue - tmp -= theta[j] * resid[l - j - 1] - - resid[l] = tmp - - if not np.isnan(tmp): - nu += 1 - ssq += tmp * tmp + return _arima.arima_transpar(params_in, arma, trans) - res = ssq / nu +# %% ../../nbs/src/arima.ipynb 19 +def arima_css(y, arma, phi, theta): + return _arima.arima_css(y, arma, phi, theta) - return res, resid - -# %% ../../nbs/src/arima.ipynb 25 -@njit(nogil=NOGIL, cache=CACHE) -def _make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(float).eps): +# %% ../../nbs/src/arima.ipynb 21 +def make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(float).eps): # check nas phi # check nas theta p = len(phi) @@ -387,14 +77,14 @@ def _make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(float).eps): if p > 0: T[:p, 0] = phi if r > 1: - for i in range(1, r): - T[i - 1, i] = 1 + idx = np.arange(1, r) + T[idx - 1, idx] = 1 if d > 0: T[r] = Z if d > 1: - for ind in range(1, d): - T[r + ind, r + ind - 1] = 1 + idx = np.arange(1, d) + T[r + idx, r + idx - 1] = 1 if q < r - 1: theta = np.concatenate((theta, np.zeros(r - 1 - q))) @@ -412,187 +102,56 @@ def _make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(float).eps): Pn[0, 0] = 1 / (1 - phi[0] ** 2) if p > 0 else 1.0 if d > 0: - for i in range(d): - Pn[r + i, r + i] = kappa - - return phi, theta, delta, Z, a, P, T, V, h, Pn + idx = np.arange(d) + Pn[r + idx, r + idx] = kappa + return { + "phi": phi, + "theta": theta, + "delta": delta, + "Z": Z, + "a": a, + "P": P, + "T": T, + "V": V, + "h": h, + "Pn": Pn, + } -def make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(np.float64).eps): - keys = ["phi", "theta", "delta", "Z", "a", "P", "T", "V", "h", "Pn"] - res = _make_arima(phi, theta, delta, kappa, tol) - return dict(zip(keys, res)) - -# %% ../../nbs/src/arima.ipynb 27 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/arima.ipynb 23 def arima_like(y, phi, theta, delta, a, P, Pn, up, use_resid): - n = len(y) - rd = len(a) - p = len(phi) - q = len(theta) - d = len(delta) - r = rd - d - - sumlog = 0.0 - ssq = 0.0 - nu = 0 - - P = P.ravel() - Pnew = Pn.ravel() - anew = np.empty(rd) - M = np.empty(rd) - if d > 0: - mm = np.empty(rd * rd) - if use_resid: - rsResid = np.empty(n) - - for l in range(n): - for i in range(r): - tmp = a[i + 1] if i < r - 1 else 0.0 - if i < p: - tmp += phi[i] * a[0] - anew[i] = tmp - if d > 0: - for i in range(r + 1, rd): - anew[i] = a[i - 1] - tmp = a[0] - for i in range(d): - tmp += delta[i] * a[r + i] - anew[r] = tmp - if l > up: - if d == 0: - for i in range(r): - vi = 0.0 - if i == 0: - vi = 1.0 - elif i - 1 < q: - vi = theta[i - 1] - for j in range(r): - tmp = 0.0 - if j == 0: - tmp = vi - elif j - 1 < q: - tmp = vi * theta[j - 1] - if i < p and j < p: - tmp += phi[i] * phi[j] * P[0] - if i < r - 1 and j < r - 1: - tmp += P[i + 1 + r * (j + 1)] - if i < p and j < r - 1: - tmp += phi[i] * P[j + 1] - if j < p and i < r - 1: - tmp += phi[j] * P[i + 1] - Pnew[i + r * j] = tmp - else: - # mm = TP - for i in range(r): - for j in range(rd): - tmp = 0.0 - if i < p: - tmp += phi[i] * P[rd * j] - if i < r - 1: - tmp += P[i + 1 + rd * j] - mm[i + rd * j] = tmp - for j in range(rd): - tmp = P[rd * j] - for k in range(d): - tmp += delta[k] * P[r + k + rd * j] - mm[r + rd * j] = tmp - for i in range(1, d): - for j in range(rd): - mm[r + i + rd * j] = P[r + i - 1 + rd * j] - - # Pnew = mmT' - for i in range(r): - for j in range(rd): - tmp = 0.0 - if i < p: - tmp += phi[i] * mm[j] - if i < r - 1: - tmp += mm[rd * (i + 1) + j] - Pnew[j + rd * i] = tmp - for j in range(rd): - tmp = mm[j] - for k in range(d): - tmp += delta[k] * mm[rd * (r + k) + j] - Pnew[rd * r + j] = tmp - for i in range(1, d): - for j in range(rd): - Pnew[rd * (r + i) + j] = mm[rd * (r + i - 1) + j] - for i in range(q + 1): - vi = 1.0 if i == 0 else theta[i - 1] - for j in range(q + 1): - Pnew[i + rd * j] += vi * (1.0 if j == 0 else theta[j - 1]) - - if not math.isnan(y[l]): - resid = y[l] - anew[0] - for i in range(d): - resid -= delta[i] * anew[r + i] - for i in range(rd): - tmp = Pnew[i] - for j in range(d): - tmp += Pnew[i + (r + j) * rd] * delta[j] - M[i] = tmp - gain = M[0] - for j in range(d): - gain += delta[j] * M[r + j] - if gain < 1e4: - nu += 1 - ssq += resid * resid / gain if gain != 0.0 else math.inf - sumlog += math.log(gain) - if use_resid: - rsResid[l] = resid / math.sqrt(gain) if gain != 0.0 else math.inf - for i in range(rd): - a[i] = anew[i] + M[i] * resid / gain if gain != 0.0 else math.inf - for i in range(rd): - for j in range(rd): - P[i + j * rd] = ( - Pnew[i + j * rd] - M[i] * M[j] / gain - if gain != 0.0 - else math.inf - ) - else: - a[:] = anew[:] - P[:] = Pnew[:] - if use_resid: - rsResid[l] = np.nan + rsResid = np.empty_like(y) + else: + rsResid = np.empty(0) + ssq, sumlog, nu = _arima.arima_like( + y, + phi, + theta, + delta, + a, + P.ravel(), + Pn.ravel(), + up, + use_resid, + rsResid, + ) if not use_resid: rsResid = None return ssq, sumlog, nu, rsResid -# %% ../../nbs/src/arima.ipynb 29 -@njit(nogil=NOGIL, cache=CACHE) -def diff1d(x, lag, differences): +# %% ../../nbs/src/arima.ipynb 25 +def diff(x, lag, differences): + x = np.asarray(x, dtype=np.float64) y = x.copy() for _ in range(differences): x = y.copy() - for i in range(lag): - y[i] = np.nan - for i in range(lag, x.size): - y[i] = x[i] - x[i - lag] - return y - - -@njit(nogil=NOGIL, cache=CACHE) -def diff2d(x, lag, differences): - y = np.empty_like(x) - for j in range(x.shape[1]): - y[:, j] = diff1d(x[:, j], lag, differences) - return y - + y[:lag] = np.nan + y[lag:] = x[lag:] - x[:-lag] + nans = lag * differences + return y[nans:] -def diff(x, lag, differences): - if x.ndim == 1: - y = diff1d(x, lag, differences) - nan_mask = np.isnan(y) - elif x.ndim == 2: - y = diff2d(x, lag, differences) - nan_mask = np.isnan(y).all(1) - else: - raise ValueError(x.ndim) - return y[~nan_mask] - -# %% ../../nbs/src/arima.ipynb 30 +# %% ../../nbs/src/arima.ipynb 26 def fixed_params_from_dict( fixed_dict: dict, order: tuple, seasonal: dict, intercept: bool, n_ex: int ): @@ -616,7 +175,7 @@ def fixed_params_from_dict( ) # prevent adding non-existing keys return list(full_dict.values()) -# %% ../../nbs/src/arima.ipynb 32 +# %% ../../nbs/src/arima.ipynb 28 def arima( x: np.ndarray, order=(0, 0, 0), @@ -634,7 +193,7 @@ def arima( optim_control={"maxiter": 100}, ): SSG = SSinit == "Gardner1980" - x = x.copy() + x = x.astype(np.float64, copy=True) def upARIMA(mod, phi, theta): p = len(phi) @@ -746,23 +305,26 @@ def maInvert(ma): # fixed # mask - arma = ( - *order[::2], - *seasonal["order"][::2], - seasonal["period"], - order[1], - seasonal["order"][1], + arma = np.array( + [ + *order[::2], + *seasonal["order"][::2], + seasonal["period"], + order[1], + seasonal["order"][1], + ], + dtype=np.intc, ) - narma = sum(arma[:4]) + narma = arma[:4].sum().item() # xtsp = init x, end x and frequency # tsp(x) = None Delta = np.array([1.0]) for i in range(order[1]): - Delta = tsconv(Delta, np.array([1.0, -1.0])) + Delta = convolve(Delta, np.array([1.0, -1.0])) for i in range(seasonal["order"][1]): - Delta = tsconv(Delta, np.array([1] + [0] * (seasonal["period"] - 1) + [-1])) + Delta = convolve(Delta, np.array([1] + [0] * (seasonal["period"] - 1) + [-1])) Delta = -Delta[1:] nd = order[1] + seasonal["order"][1] n_used = (~np.isnan(x)).sum() - len(Delta) @@ -889,7 +451,7 @@ def arma_css_op(p, x): if ncxreg > 0: x -= np.dot(xreg, par[narma + np.arange(ncxreg)]) - res, resid = arima_css(x, arma, phi, theta, ncond) + res, _ = arima_css(x, arma, phi, theta) if math.isinf(res): import sys @@ -922,7 +484,7 @@ def arma_css_op(p, x): mod = make_arima(phi, theta, Delta, kappa) if ncxreg > 0: x -= np.dot(xreg, coef[narma + np.arange(ncxreg)]) - val = arima_css(x, arma, phi, theta, ncond) + val = arima_css(x, arma, phi, theta) sigma2 = val[0] var = None if no_optim else res.hess_inv / n_used else: @@ -939,7 +501,7 @@ def arma_css_op(p, x): # only update the initial parameters if they're valid candidate = init.copy() candidate[mask] = res.x - phi, _ = arima_transpar(candidate, arma, False) + phi, _ = arima_transpar(candidate, arma, transform_pars) if np.logical_and(phi > -math.pi / 2, phi < math.pi / 2).all(): init = candidate if arma[0] > 0: @@ -950,13 +512,16 @@ def arma_css_op(p, x): raise ValueError("non-stationary seasonal AR part from CSS") ncond = 0 if transform_pars: - init = ARIMA_invtrans(init, arma) - if arma[1] > 0: - ind = arma[0] + np.arange(arma[1]) - init[ind] = maInvert(init[ind]) - if arma[3] > 0: - ind = np.sum(arma[:3]) + np.arange(arma[3]) - init[ind] = maInvert(init[ind]) + candidate = ARIMA_invtrans(init, arma) + phi, _ = arima_transpar(candidate, arma, transform_pars) + if np.logical_and(phi > -math.pi / 2, phi < math.pi / 2).all(): + init = candidate + if arma[1] > 0: + ind = arma[0] + np.arange(arma[1]) + init[ind] = maInvert(init[ind]) + if arma[3] > 0: + ind = np.sum(arma[:3]) + np.arange(arma[3]) + init[ind] = maInvert(init[ind]) trarma = arima_transpar(init, arma, transform_pars) mod = make_arima(trarma[0], trarma[1], Delta, kappa, SSinit) if no_optim: @@ -1048,7 +613,7 @@ def arma_css_op(p, x): "mask": mask, "loglik": -0.5 * value, "aic": aic, - "arma": arma, + "arma": tuple(x.item() for x in arma), "residuals": resid, #'series': series, "code": res.status, @@ -1058,55 +623,28 @@ def arma_css_op(p, x): } return ans -# %% ../../nbs/src/arima.ipynb 40 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/arima.ipynb 36 def kalman_forecast(n, Z, a, P, T, V, h): - p = len(a) - a = a.copy() - anew = np.empty(p) - Pnew = np.empty((p, p)) - mm = np.empty((p, p)) + P = P.copy() forecasts = np.empty(n) se = np.empty(n) - P = P.copy() - + z = Z.reshape(-1, 1) * Z.reshape(1, -1) for l in range(n): - anew = T @ a - - a[:] = anew[:] - forecasts[l] = anew @ Z - - for i in range(p): - for j in range(p): - tmp = 0.0 - for k in range(p): - tmp += T[i, k] * P[k, j] - mm[i, j] = tmp - - for i in range(p): - for j in range(p): - tmp = V[i, j] - for k in range(p): - tmp += mm[i, k] * T[j, k] - Pnew[i, j] = tmp - - tmp = h - for i in range(p): - for j in range(p): - P[i, j] = Pnew[i, j] - tmp += Z[i] * Z[j] * P[i, j] - se[l] = tmp - + a = T @ a + forecasts[l] = a @ Z + mm = T @ P + P = V + mm @ T.T + se[l] = h + np.sum(z * P) return forecasts, se -# %% ../../nbs/src/arima.ipynb 43 +# %% ../../nbs/src/arima.ipynb 39 def checkarima(obj): if obj["var_coef"] is None: return False return any(np.isnan(np.sqrt(np.diag(obj["var_coef"])))) -# %% ../../nbs/src/arima.ipynb 44 +# %% ../../nbs/src/arima.ipynb 40 def predict_arima(model, n_ahead, newxreg=None, se_fit=True): myNCOL = lambda x: x.shape[1] if x is not None else 0 @@ -1165,7 +703,7 @@ def predict_arima(model, n_ahead, newxreg=None, se_fit=True): return pred -# %% ../../nbs/src/arima.ipynb 48 +# %% ../../nbs/src/arima.ipynb 44 def convert_coef_name(name, inverse=False): if not inverse: if "ex" in name: @@ -1187,13 +725,13 @@ def convert_coef_name(name, inverse=False): else: return name -# %% ../../nbs/src/arima.ipynb 49 +# %% ../../nbs/src/arima.ipynb 45 def change_drift_name(model_coef, inverse=False): return { convert_coef_name(name, inverse): value for name, value in model_coef.items() } -# %% ../../nbs/src/arima.ipynb 50 +# %% ../../nbs/src/arima.ipynb 46 def myarima( x, order=(0, 0, 0), @@ -1292,7 +830,7 @@ def myarima( raise e return {"ic": math.inf} -# %% ../../nbs/src/arima.ipynb 53 +# %% ../../nbs/src/arima.ipynb 49 def search_arima( x, d=0, @@ -1386,7 +924,7 @@ def search_arima( ) return best_fit -# %% ../../nbs/src/arima.ipynb 55 +# %% ../../nbs/src/arima.ipynb 51 def arima2(x, model, xreg, method): m = model["arma"][4] # 5 use_drift = "drift" in model["coef"].keys() @@ -1453,7 +991,7 @@ def arima2(x, model, xreg, method): refit["coef"] = change_drift_name(refit["coef"]) return refit -# %% ../../nbs/src/arima.ipynb 56 +# %% ../../nbs/src/arima.ipynb 52 def Arima( x, order=(0, 0, 0), @@ -1547,7 +1085,7 @@ def Arima( tmp["sigma2"] = np.nansum(tmp["residuals"] ** 2) / (nstar - npar + 1) return tmp -# %% ../../nbs/src/arima.ipynb 64 +# %% ../../nbs/src/arima.ipynb 60 def arima_string(model, padding=False): order = tuple(model["arma"][i] for i in [0, 5, 1, 2, 6, 3, 4]) m = order[6] @@ -1577,11 +1115,11 @@ def arima_string(model, padding=False): return result -# %% ../../nbs/src/arima.ipynb 67 +# %% ../../nbs/src/arima.ipynb 63 def is_constant(x): return np.all(x[0] == x) -# %% ../../nbs/src/arima.ipynb 68 +# %% ../../nbs/src/arima.ipynb 64 def forecast_arima( model, h=None, @@ -1676,7 +1214,7 @@ def forecast_arima( return ans -# %% ../../nbs/src/arima.ipynb 75 +# %% ../../nbs/src/arima.ipynb 71 def fitted_arima(model, h=1): """Returns h-step forecasts for the data used in fitting the model.""" if h == 1: @@ -1692,7 +1230,7 @@ def fitted_arima(model, h=1): else: raise NotImplementedError("h > 1") -# %% ../../nbs/src/arima.ipynb 80 +# %% ../../nbs/src/arima.ipynb 76 def seas_heuristic(x, period): # nperiods = period > 1 season = math.nan @@ -1704,7 +1242,7 @@ def seas_heuristic(x, period): season = max(0, min(1, 1 - vare / np.var(remainder + seasonal, ddof=1))) return season -# %% ../../nbs/src/arima.ipynb 82 +# %% ../../nbs/src/arima.ipynb 78 def nsdiffs(x, test="seas", alpha=0.05, period=1, max_D=1, **kwargs): D = 0 if alpha < 0.01: @@ -1767,7 +1305,7 @@ def run_tests(x, test, alpha): dodiff = False return D -# %% ../../nbs/src/arima.ipynb 84 +# %% ../../nbs/src/arima.ipynb 80 def ndiffs(x, alpha=0.05, test="kpss", kind="level", max_d=2): x = x[~np.isnan(x)] d = 0 @@ -1812,13 +1350,13 @@ def run_tests(x, test, alpha): return d - 1 return d -# %% ../../nbs/src/arima.ipynb 86 +# %% ../../nbs/src/arima.ipynb 82 def newmodel(p, d, q, P, D, Q, constant, results): curr = np.array([p, d, q, P, D, Q, constant]) in_results = (curr == results[:, :7]).all(1).any() return not in_results -# %% ../../nbs/src/arima.ipynb 88 +# %% ../../nbs/src/arima.ipynb 84 def auto_arima_f( x, d=None, @@ -2356,11 +1894,11 @@ def try_params(p, d, q, P, D, Q, constant, k, bestfit): return bestfit -# %% ../../nbs/src/arima.ipynb 90 +# %% ../../nbs/src/arima.ipynb 86 def forward_arima(fitted_model, y, xreg=None, method="CSS-ML"): return Arima(x=y, model=fitted_model, xreg=xreg, method=method) -# %% ../../nbs/src/arima.ipynb 99 +# %% ../../nbs/src/arima.ipynb 95 def print_statsforecast_ARIMA(model, digits=3, se=True): print(arima_string(model, padding=False)) if model["lambda"] is not None: @@ -2390,7 +1928,7 @@ def print_statsforecast_ARIMA(model, digits=3, se=True): if not np.isnan(model["aic"]): print(f'AIC={round(model["aic"], 2)}') -# %% ../../nbs/src/arima.ipynb 101 +# %% ../../nbs/src/arima.ipynb 97 class ARIMASummary: """ARIMA Summary.""" @@ -2403,7 +1941,7 @@ def __repr__(self): def summary(self): return print_statsforecast_ARIMA(self.model) -# %% ../../nbs/src/arima.ipynb 102 +# %% ../../nbs/src/arima.ipynb 98 class AutoARIMA: """An AutoARIMA estimator. diff --git a/setup.py b/setup.py index 34987d13b..4bd43e74a 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ import setuptools from configparser import ConfigParser -from pybind11.setup_helpers import Pybind11Extension +from pybind11.setup_helpers import ParallelCompile, Pybind11Extension, naive_recompile # note: all settings are in settings.ini; edit there, not here config = ConfigParser(delimiters=['=']) @@ -61,6 +61,7 @@ cxx_std=17, ) ] +ParallelCompile("CMAKE_BUILD_PARALLEL_LEVEL", needs_recompile=naive_recompile).install() setuptools.setup( name = 'statsforecast', diff --git a/src/arima.cpp b/src/arima.cpp new file mode 100644 index 000000000..4d3879a49 --- /dev/null +++ b/src/arima.cpp @@ -0,0 +1,584 @@ +#include +#include +#include +#include + +#include +#include + +namespace arima { +namespace py = pybind11; + +void partrans(int p, const double *raw, double *newv) { + std::transform(raw, raw + p, newv, [](double x) { return std::tanh(x); }); + std::vector work(newv, newv + p); + for (int j = 1; j < p; ++j) { + for (int k = 0; k < j; ++k) { + work[k] -= newv[j] * newv[j - k - 1]; + } + std::copy(work.begin(), work.begin() + j, newv); + } +} + +std::tuple, py::array_t> +arima_transpar(const py::array_t params_inv, + const py::array_t armav, bool trans) { + auto arma = armav.data(); + auto params_in = params_inv.data(); + int mp = arma[0]; + int mq = arma[1]; + int msp = arma[2]; + int msq = arma[3]; + int ns = arma[4]; + int p = mp + ns * msp; + int q = mq + ns * msq; + auto params = std::vector(params_in, params_in + params_inv.size()); + py::array_t phiv(p); + py::array_t thetav(q); + auto phi = phiv.mutable_data(); + auto theta = thetav.mutable_data(); + if (trans) { + if (mp > 0) { + partrans(mp, params_in, params.data()); + } + int v = mp + mq; + if (msp > 0) { + partrans(msp, params_in + v, params.data() + v); + } + } + if (ns > 0) { + std::copy(params.begin(), params.begin() + mp, phi); + std::fill(phi + mp, phi + p, 0.0); + std::copy(params.begin() + mp, params.begin() + mp + mq, theta); + std::fill(theta + mq, theta + q, 0.0); + for (int j = 0; j < msp; ++j) { + phi[(j + 1) * ns - 1] += params[j + mp + mq]; + for (int i = 0; i < mp; ++i) { + phi[(j + 1) * ns + i] -= params[i] * params[j + mp + mq]; + } + } + for (int j = 0; j < msq; ++j) { + theta[(j + 1) * ns - 1] += params[j + mp + mq + msp]; + for (int i = 0; i < mq; ++i) { + theta[(j + 1) * ns + i] += params[i + mp] * params[j + mp + mq + msp]; + } + } + } else { + std::copy(params.begin(), params.begin() + mp, phi); + std::copy(params.begin() + mp, params.begin() + mp + mq, theta); + } + return {phiv, thetav}; +} + +std::tuple> +arima_css(const py::array_t yv, const py::array_t armav, + const py::array_t phiv, const py::array_t thetav) { + int n = static_cast(yv.size()); + int p = static_cast(phiv.size()); + int q = static_cast(thetav.size()); + auto y = yv.data(); + auto arma = armav.data(); + auto phi = phiv.data(); + auto theta = thetav.data(); + int ncond = arma[0] + arma[5] + arma[4] * (arma[2] + arma[6]); + int nu = 0; + double ssq = 0.0; + + auto residv = py::array_t(n); + auto resid = residv.mutable_data(); + auto w = std::vector(y, y + yv.size()); + for (int _ = 0; _ < arma[5]; ++_) { + for (int l = n - 1; l > 0; --l) { + w[l] -= w[l - 1]; + } + } + int ns = arma[4]; + for (int _ = 0; _ < arma[6]; ++_) { + for (int l = n - 1; l >= ns; --l) { + w[l] -= w[l - ns]; + } + } + for (int l = ncond; l < n; ++l) { + double tmp = w[l]; + for (int j = 0; j < p; ++j) { + tmp -= phi[j] * w[l - j - 1]; + } + for (int j = 0; j < std::min(l - ncond, q); ++j) { + if (l - j - 1 < 0) { + continue; + } + tmp -= theta[j] * resid[l - j - 1]; + } + resid[l] = tmp; + if (!std::isnan(tmp)) { + nu++; + ssq += tmp * tmp; + } + } + return {ssq / nu, residv}; +} + +std::tuple +arima_like(const py::array_t yv, const py::array_t phiv, + const py::array_t thetav, const py::array_t deltav, + py::array_t av, py::array_t Pv, + py::array_t Pnewv, int up, bool use_resid, + py::array_t rsResidv) { + int n = static_cast(yv.size()); + int d = static_cast(deltav.size()); + int rd = static_cast(av.size()); + int p = static_cast(phiv.size()); + int q = static_cast(thetav.size()); + auto y = yv.data(); + auto phi = phiv.data(); + auto theta = thetav.data(); + auto delta = deltav.data(); + auto a = av.mutable_data(); + auto P = Pv.mutable_data(); + auto Pnew = Pnewv.mutable_data(); + auto rsResid = rsResidv.mutable_data(); + double ssq = 0.0; + double sumlog = 0.0; + int nu = 0; + int r = rd - d; + + std::vector anew(rd); + std::vector M(rd); + std::vector mm; + if (d > 0) { + mm.resize(rd * rd); + } + double tmp; + for (int l = 0; l < n; ++l) { + for (int i = 0; i < r; ++i) { + if (i < r - 1) { + tmp = a[i + 1]; + } else { + tmp = 0.0; + } + if (i < p) { + tmp += phi[i] * a[0]; + } + anew[i] = tmp; + } + if (d > 0) { + for (int i = r + 1; i < rd; ++i) { + anew[i] = a[i - 1]; + } + tmp = a[0]; + for (int i = 0; i < d; ++i) { + tmp += delta[i] * a[r + i]; + } + anew[r] = tmp; + } + if (l > up) { + if (d == 0) { + for (int i = 0; i < r; ++i) { + double vi = 0.0; + if (i == 0) { + vi = 1.0; + } else if (i - 1 < q) { + vi = theta[i - 1]; + } + for (int j = 0; j < r; ++j) { + tmp = 0.0; + if (j == 0) { + tmp = vi; + } else if (j - 1 < q) { + tmp = vi * theta[j - 1]; + } + if (i < p && j < p) { + tmp += phi[i] * phi[j] * P[0]; + } + if (i < r - 1 && j < r - 1) { + tmp += P[i + 1 + r * (j + 1)]; + } + if (i < p && j < r - 1) { + tmp += phi[i] * P[j + 1]; + } + if (j < p && i < r - 1) { + tmp += phi[j] * P[i + 1]; + } + Pnew[i + r * j] = tmp; + } + } + } else { + for (int i = 0; i < r; ++i) { + for (int j = 0; j < rd; ++j) { + tmp = 0.0; + if (i < p) { + tmp += phi[i] * P[rd * j]; + } + if (i < r - 1) { + tmp += P[i + 1 + rd * j]; + } + mm[i + rd * j] = tmp; + } + } + for (int j = 0; j < rd; ++j) { + tmp = P[rd * j]; + for (int k = 0; k < d; ++k) { + tmp += delta[k] * P[r + k + rd * j]; + } + mm[r + rd * j] = tmp; + } + for (int i = 1; i < d; ++i) { + for (int j = 0; j < rd; ++j) { + mm[r + i + rd * j] = P[r + i - 1 + rd * j]; + } + } + for (int i = 0; i < r; ++i) { + for (int j = 0; j < rd; ++j) { + tmp = 0.0; + if (i < p) { + tmp += phi[i] * mm[j]; + } + if (i < r - 1) { + tmp += mm[rd * (i + 1) + j]; + } + Pnew[j + rd * i] = tmp; + } + } + for (int j = 0; j < rd; ++j) { + tmp = mm[j]; + for (int k = 0; k < d; ++k) { + tmp += delta[k] * mm[rd * (r + k) + j]; + } + Pnew[rd * r + j] = tmp; + } + for (int i = 1; i < d; ++i) { + for (int j = 0; j < rd; ++j) { + Pnew[rd * (r + i) + j] = mm[rd * (r + i - 1) + j]; + } + } + for (int i = 0; i < q + 1; ++i) { + double vi; + if (i == 0) { + vi = 1.0; + } else { + vi = theta[i - 1]; + } + for (int j = 0; j < q + 1; ++j) { + if (j == 0) { + Pnew[i + rd * j] += vi; + } else { + Pnew[i + rd * j] += vi * theta[j - 1]; + } + } + } + } + } + if (!std::isnan(y[l])) { + double resid = y[l] - anew[0]; + for (int i = 0; i < d; ++i) { + resid -= delta[i] * anew[r + i]; + } + for (int i = 0; i < rd; ++i) { + tmp = Pnew[i]; + for (int j = 0; j < d; ++j) { + tmp += Pnew[i + (r + j) * rd] * delta[j]; + } + M[i] = tmp; + } + double gain = M[0]; + for (int j = 0; j < d; ++j) { + gain += delta[j] * M[r + j]; + } + if (gain < 1e4) { + nu++; + if (gain == 0) { + ssq = std::numeric_limits::infinity(); + } else { + ssq += resid * resid / gain; + } + sumlog += std::log(gain); + } + if (use_resid) { + if (gain == 0) { + rsResid[l] = std::numeric_limits::infinity(); + } else { + rsResid[l] = resid / std::sqrt(gain); + } + } + if (gain == 0) { + for (int i = 0; i < rd; ++i) { + a[i] = std::numeric_limits::infinity(); + for (int j = 0; j < rd; ++j) { + Pnew[i + j * rd] = std::numeric_limits::infinity(); + } + } + } else { + for (int i = 0; i < rd; ++i) { + a[i] = anew[i] + M[i] * resid / gain; + for (int j = 0; j < rd; ++j) { + P[i + j * rd] = Pnew[i + j * rd] - M[i] * M[j] / gain; + } + } + } + } else { + std::copy(anew.begin(), anew.end(), a); + std::copy(Pnew, Pnew + rd * rd, P); + if (use_resid) { + rsResid[l] = std::numeric_limits::quiet_NaN(); + } + } + } + return {ssq, sumlog, nu}; +} + +void inclu2(int np, const double *xnext, double *xrow, double ynext, double *d, + double *rbar, double *thetab) { + std::copy(xnext, xnext + np, xrow); + int ithisr = 0; + for (int i = 0; i < np; ++i) { + if (xrow[i] != 0.0) { + double xi = xrow[i]; + double di = d[i]; + double dpi = di + xi * xi; + d[i] = dpi; + double cbar, sbar; + if (dpi == 0) { + cbar = std::numeric_limits::infinity(); + sbar = std::numeric_limits::infinity(); + } else { + cbar = di / dpi; + sbar = xi / dpi; + } + for (int k = i + 1; k < np; ++k) { + double xk = xrow[k]; + double rbthis = rbar[ithisr]; + xrow[k] = xk - xi * rbthis; + rbar[ithisr++] = cbar * rbthis + sbar * xk; + } + double xk = ynext; + ynext = xk - xi * thetab[i]; + thetab[i] = cbar * thetab[i] + sbar * xk; + if (di == 0.0) { + return; + } + } else { + ithisr += np - i - 1; + } + } +} + +void getQ0(const py::array_t phiv, const py::array_t thetav, + py::array_t resv) { + auto phi = phiv.data(); + auto theta = thetav.data(); + auto res = resv.mutable_data(); + int p = static_cast(phiv.size()); + int q = static_cast(thetav.size()); + int r = std::max(p, q + 1); + int np = r * (r + 1) / 2; + int nrbar = np * (np - 1) / 2; + int ind = 0; + + std::vector V(np); + for (int j = 0; j < r; ++j) { + double vj = 0.0; + if (j == 0) { + vj = 1.0; + } else if (j - 1 < q) { + vj = theta[j - 1]; + } + for (int i = j; i < r; ++i) { + double vi = 0.0; + if (i == 0) { + vi = 1.0; + } else if (i - 1 < q) { + vi = theta[i - 1]; + } + V[ind++] = vi * vj; + } + } + if (r == 1) { + if (p == 0) { + res[0] = 1.0; + } else { + res[0] = 1.0 / (1 - phi[0] * phi[0]); + } + return; + } + if (p > 0) { + std::vector rbar(nrbar); + std::vector thetab(np); + std::vector xnext(np); + std::vector xrow(np); + ind = 0; + int ind1 = -1; + int npr = np - r; + int npr1 = npr + 1; + int indj = npr; + int ind2 = npr - 1; + for (int j = 0; j < r; ++j) { + double phij = j < p ? phi[j] : 0.0; + xnext[indj++] = 0.0; + int indi = npr1 + j; + for (int i = j; i < r; ++i) { + double ynext = V[ind++]; + double phii = i < p ? phi[i] : 0.0; + if (j != r - 1) { + xnext[indj] = -phii; + if (i != r - 1) { + xnext[indi] -= phij; + xnext[++ind1] = -1.0; + } + } + xnext[npr] = -phii * phij; + if (++ind2 >= np) { + ind2 = 0; + } + xnext[ind2] += 1.0; + inclu2(np, xnext.data(), xrow.data(), ynext, res, rbar.data(), + thetab.data()); + xnext[ind2] = 0.0; + if (i != r - 1) { + xnext[indi++] = 0.0; + xnext[ind1] = 0.0; + } + } + } + int ithisr = nrbar - 1; + int im = np - 1; + for (int i = 0; i < np; ++i) { + double bi = thetab[im]; + int jm = np - 1; + for (int j = 0; j < i; ++j) { + bi -= rbar[ithisr--] * res[jm--]; + } + res[im--] = bi; + } + ind = npr; + for (int i = 0; i < r; ++i) { + xnext[i] = res[ind++]; + } + ind = np - 1; + ind1 = npr - 1; + for (int i = 0; i < npr; ++i) { + res[ind--] = res[ind1--]; + } + std::copy(xnext.begin(), xnext.begin() + r, res); + } else { + int indn = np; + ind = np; + for (int i = 0; i < r; ++i) { + for (int j = 0; j < i + 1; ++j) { + --ind; + res[ind] = V[ind]; + if (j != 0) { + res[ind] += res[--indn]; + } + } + } + } + ind = np; + for (int i = r - 1; i > 0; --i) { + for (int j = r - 1; j > i - 1; --j) { + res[r * i + j] = res[--ind]; + } + } + for (int i = 0; i < r - 1; ++i) { + for (int j = i + 1; j < r; ++j) { + res[i + r * j] = res[j + r * i]; + } + } +} + +py::array_t arima_gradtrans(const py::array_t xv, + const py::array_t armav) { + constexpr double eps = 1e-3; + auto x = xv.data(); + auto arma = armav.data(); + int n = static_cast(xv.size()); + int mp = arma[0]; + int mq = arma[1]; + int msp = arma[2]; + + auto w1 = std::array(); + auto w2 = std::array(); + auto w3 = std::array(); + py::array_t outv({n, n}); + auto out = outv.mutable_data(); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + out[i * n + j] = (i == j) ? 1.0 : 0.0; + } + } + if (mp > 0) { + std::copy(x, x + mp, w1.begin()); + partrans(mp, w1.data(), w2.data()); + for (int i = 0; i < mp; ++i) { + w1[i] += eps; + partrans(mp, w1.data(), w3.data()); + for (int j = 0; j < mp; ++j) { + out[n * i + j] = (w3[j] - w2[j]) / eps; + } + w1[i] -= eps; + } + } + if (msp > 0) { + int v = mp + mq; + std::copy(x + v, x + v + msp, w1.begin()); + partrans(msp, w1.data(), w2.data()); + for (int i = 0; i < msp; ++i) { + w1[i] += eps; + partrans(msp, w1.data(), w3.data()); + for (int j = 0; j < msp; ++j) { + out[n * (i + v) + j + v] = (w3[j] - w2[j]) / eps; + } + w1[i] -= eps; + } + } + return outv; +} + +py::array_t arima_undopars(const py::array_t xv, + const py::array_t armav) { + auto x = xv.data(); + auto arma = armav.data(); + int mp = arma[0]; + int mq = arma[1]; + int msp = arma[2]; + py::array_t outv{xv.size()}; + auto out = outv.mutable_data(); + std::copy(xv.data(), xv.data() + xv.size(), out); + if (mp > 0) { + partrans(mp, x, out); + } + int v = mp + mq; + if (msp > 0) { + partrans(msp, x + v, out + v); + } + return outv; +} + +void invpartrans(int p, const py::array_t phiv, + py::array_t outv) { + auto phi = phiv.data(); + auto out = outv.mutable_data(); + std::copy(phi, phi + p, out); + std::vector work(phi, phi + p); + for (int j = p - 1; j > 0; --j) { + double a = out[j]; + for (int k = 0; k < j; ++k) { + work[k] = (out[k] + a * out[j - k - 1]) / (1 - a * a); + } + std::copy(work.begin(), work.begin() + j, out); + } + for (int j = 0; j < p; ++j) { + out[j] = std::atanh(out[j]); + } +} + +void init(py::module_ &m) { + py::module_ arima = m.def_submodule("arima"); + arima.def("arima_css", &arima_css); + arima.def("arima_like", &arima_like); + arima.def("getQ0", &getQ0); + arima.def("arima_gradtrans", &arima_gradtrans); + arima.def("arima_undopars", &arima_undopars); + arima.def("invpartrans", &invpartrans); + arima.def("arima_transpar", &arima_transpar); +} +} // namespace arima diff --git a/src/statsforecast.cpp b/src/statsforecast.cpp index f240b4a1f..edbbed6b5 100644 --- a/src/statsforecast.cpp +++ b/src/statsforecast.cpp @@ -2,10 +2,18 @@ namespace py = pybind11; -namespace ets { -void init(py::module_ &); +namespace ets +{ + void init(py::module_ &); } -PYBIND11_MODULE(_lib, m) { +namespace arima +{ + void init(py::module_ &); +} + +PYBIND11_MODULE(_lib, m) +{ + arima::init(m); ets::init(m); }