fractional_differencing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "89a52074-aa13-458f-bc22-e3d76c301b0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "id": "f7ea0850-6d19-419e-b671-13400664cb61",
   "metadata": {},
   "outputs": [],
   "source": [
    "def frac_diff(x, d):\n",
    "    \"\"\"\n",
    "    Fractionally difference time series\n",
    "\n",
    "    :param x: numeric vector or univariate time series\n",
    "    :param d: number specifying the fractional difference order.\n",
    "    :return: fractionally differenced series\n",
    "    \"\"\"\n",
    "    if np.isnan(np.sum(x)):\n",
    "        return None\n",
    "\n",
    "    n = len(x)\n",
    "    if n < 2:\n",
    "        return None\n",
    "\n",
    "    x = np.subtract(x, np.mean(x))\n",
    "\n",
    "    # calculate weights\n",
    "    weights = [0] * n\n",
    "    weights[0] = -d\n",
    "    for k in range(2, n):\n",
    "        weights[k - 1] = weights[k - 2] * (k - 1 - d) / k\n",
    "\n",
    "    # difference series\n",
    "    ydiff = list(x)\n",
    "\n",
    "    for i in range(0, n):\n",
    "        dat = x[:i]\n",
    "        w = weights[:i]\n",
    "        ydiff[i] = x[i] + np.dot(w, dat[::-1])\n",
    "\n",
    "    return ydiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "a65c7ad5-5b97-42f2-a9ce-73b6c4785d22",
   "metadata": {},
   "outputs": [],
   "source": [
    "def frac_diff_vect(x, d):\n",
    "    \"\"\"\n",
    "    Fractionally difference time series\n",
    "\n",
    "    :param x: numeric vector or univariate time series\n",
    "    :param d: number specifying the fractional difference order\n",
    "    :return: fractionally differenced series\n",
    "    \"\"\"\n",
    "    if np.isnan(np.sum(x)):\n",
    "        return None\n",
    "\n",
    "    n = len(x)\n",
    "    if n < 2:\n",
    "        return None\n",
    "\n",
    "    x = np.subtract(x, np.mean(x))\n",
    "\n",
    "    # calculate weights\n",
    "    num = -d*np.cumprod(np.arange(1-d, n-1-d, 1))\n",
    "    den = np.cumprod(np.arange(2,n,1))\n",
    "    weights = np.divide(num, den)\n",
    "    weights = np.insert(weights, 0, -d)\n",
    "    weights = np.insert(weights, weights.shape, 0)\n",
    "\n",
    "    # difference series\n",
    "    ydiff = x.reshape(1,-1) * weights.reshape(-1,1)\n",
    "    \n",
    "    def tracify(k): \n",
    "        return np.trace(np.fliplr(ydiff), k)\n",
    "    \n",
    "    vfunc = np.vectorize(tracify)\n",
    "    ydiff = x + vfunc(np.arange(x.shape[0], 0, -1))\n",
    "\n",
    "    return ydiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "id": "e06276aa-4db6-45c1-9f8f-de3be67be83f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -6.14285714,  -2.07142857,  -0.80357143,  -0.04464286,\n",
       "       -11.4921875 ,  15.94921875,  23.32128906])"
      ]
     },
     "execution_count": 190,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "frac_diff_vect(np.array([1,2,3,4,-7,15, 32]), 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "id": "bd524052-28d4-4a63-bfb1-996ef3782be6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-6.142857142857143,\n",
       " -2.0714285714285716,\n",
       " -0.8035714285714288,\n",
       " -0.044642857142857206,\n",
       " -11.4921875,\n",
       " 15.94921875,\n",
       " 23.3212890625]"
      ]
     },
     "execution_count": 191,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "frac_diff(np.array([1,2,3,4,-7,15, 32]), 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "id": "7ba398de-2cb4-40c2-ac17-22ea449a7628",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.rand(30000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "id": "f84b4181-a010-4656-8424-2a778a4031be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "34.66284656524658\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "y = frac_diff(np.array(a), 0.05)\n",
    "end = time.time()\n",
    "print(end - start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "2439ddb0-c6e8-4d76-82b4-7bdcd1998d14",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\mmariantoni\\AppData\\Local\\Temp\\1\\ipykernel_27280\\1325045861.py:21: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  weights = np.divide(num, den)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.000987768173218\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "x = frac_diff_vect(np.array(a), 0.05)\n",
    "end = time.time()\n",
    "print(end - start)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}