DM_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dm_test(actual_lst, pred1_lst, pred2_lst, h=1, crit=\"MSE\", power=2):\n",
    "    \"\"\"\n",
    "    这个是DM检验的代码\n",
    "    :param actual_lst: 真实的序列值\n",
    "    :param pred1_lst: 第一个模型预测的结果\n",
    "    :param pred2_lst: 第二个模型预测的结果\n",
    "    :param h: 预测模型是几步预测，h就是几\n",
    "    :param crit: 计算连个模型的预测偏差，的差值 d 时，使用的公式：有MSE,MAD,MAPE,poly，推荐MSE\n",
    "    :param power: 只有crit=poly是用这个，计算d时使用： (模型1的偏差)的power次方 - (模型2的偏差)的power次方\n",
    "    :return:\n",
    "    \"\"\"\n",
    "\n",
    "    # Routine for checking errors\n",
    "    def error_check():\n",
    "        rt = 0\n",
    "        msg = \"\"\n",
    "        # Check if h is an integer\n",
    "        if not isinstance(h, int):\n",
    "            rt = -1\n",
    "            msg = \"The type of the number of steps ahead (h) is not an integer.\"\n",
    "            return rt, msg\n",
    "        # Check the range of h\n",
    "        if h < 1:\n",
    "            rt = -1\n",
    "            msg = \"The number of steps ahead (h) is not large enough.\"\n",
    "            return rt, msg\n",
    "        len_act = len(actual_lst)\n",
    "        len_p1 = len(pred1_lst)\n",
    "        len_p2 = len(pred2_lst)\n",
    "        # Check if lengths of actual values and predicted values are equal\n",
    "        if len_act != len_p1 or len_p1 != len_p2 or len_act != len_p2:\n",
    "            rt = -1\n",
    "            msg = \"Lengths of actual_lst, pred1_lst and pred2_lst do not match.\"\n",
    "            return rt, msg\n",
    "        # Check range of h\n",
    "        if h >= len_act:\n",
    "            rt = -1\n",
    "            msg = \"The number of steps ahead is too large.\"\n",
    "            return rt, msg\n",
    "        # Check if criterion supported\n",
    "        if crit != \"MSE\" and crit != \"MAPE\" and crit != \"MAD\" and crit != \"poly\":\n",
    "            rt = -1\n",
    "            msg = \"The criterion is not supported.\"\n",
    "            return rt, msg\n",
    "            # Check if every value of the input lists are numerical values\n",
    "        from re import compile as re_compile\n",
    "        comp = re_compile(\"^\\d+?\\.\\d+?$\")\n",
    "\n",
    "        def compiled_regex(s):\n",
    "            \"\"\" Returns True is string is a number. \"\"\"\n",
    "            if comp.match(s) is None:\n",
    "                return s.isdigit()\n",
    "            return True\n",
    "\n",
    "        for actual, pred1, pred2 in zip(actual_lst, pred1_lst, pred2_lst):\n",
    "            is_actual_ok = compiled_regex(str(abs(actual)))\n",
    "            is_pred1_ok = compiled_regex(str(abs(pred1)))\n",
    "            is_pred2_ok = compiled_regex(str(abs(pred2)))\n",
    "            if not (is_actual_ok and is_pred1_ok and is_pred2_ok):\n",
    "                msg = \"An element in the actual_lst, pred1_lst or pred2_lst is not numeric.\"\n",
    "                rt = -1\n",
    "                return rt, msg\n",
    "        return rt, msg\n",
    "\n",
    "    # Error check\n",
    "    error_code = error_check()\n",
    "    # Raise error if cannot pass error check\n",
    "    if error_code[0] == -1:\n",
    "        raise SyntaxError(error_code[1])\n",
    "        return\n",
    "    # Import libraries\n",
    "    from scipy.stats import t\n",
    "    import collections\n",
    "    import pandas as pd\n",
    "    import numpy as np\n",
    "\n",
    "    # Initialise lists\n",
    "    e1_lst = []\n",
    "    e2_lst = []\n",
    "    d_lst = []\n",
    "\n",
    "    # convert every value of the lists into real values\n",
    "    actual_lst = pd.Series(actual_lst).apply(lambda x: float(x)).tolist()\n",
    "    pred1_lst = pd.Series(pred1_lst).apply(lambda x: float(x)).tolist()\n",
    "    pred2_lst = pd.Series(pred2_lst).apply(lambda x: float(x)).tolist()\n",
    "\n",
    "    # Length of lists (as real numbers)\n",
    "    T = float(len(actual_lst))\n",
    "\n",
    "    # construct d according to crit\n",
    "    if crit == \"MSE\":\n",
    "        for actual, p1, p2 in zip(actual_lst, pred1_lst, pred2_lst):\n",
    "            e1_lst.append((actual - p1) ** 2)\n",
    "            e2_lst.append((actual - p2) ** 2)\n",
    "        for e1, e2 in zip(e1_lst, e2_lst):\n",
    "            d_lst.append(e1 - e2)\n",
    "    elif crit == \"MAD\":\n",
    "        for actual, p1, p2 in zip(actual_lst, pred1_lst, pred2_lst):\n",
    "            e1_lst.append(abs(actual - p1))\n",
    "            e2_lst.append(abs(actual - p2))\n",
    "        for e1, e2 in zip(e1_lst, e2_lst):\n",
    "            d_lst.append(e1 - e2)\n",
    "    elif crit == \"MAPE\":\n",
    "        for actual, p1, p2 in zip(actual_lst, pred1_lst, pred2_lst):\n",
    "            e1_lst.append(abs((actual - p1) / actual))\n",
    "            e2_lst.append(abs((actual - p2) / actual))\n",
    "        for e1, e2 in zip(e1_lst, e2_lst):\n",
    "            d_lst.append(e1 - e2)\n",
    "    elif crit == \"poly\":\n",
    "        for actual, p1, p2 in zip(actual_lst, pred1_lst, pred2_lst):\n",
    "            e1_lst.append(((actual - p1)) ** (power))\n",
    "            e2_lst.append(((actual - p2)) ** (power))\n",
    "        for e1, e2 in zip(e1_lst, e2_lst):\n",
    "            d_lst.append(e1 - e2)\n",
    "\n",
    "            # Mean of d\n",
    "    mean_d = pd.Series(d_lst).mean()\n",
    "\n",
    "    # Find autocovariance and construct DM test statistics\n",
    "    def autocovariance(Xi, N, k, Xs):\n",
    "        autoCov = 0\n",
    "        T = float(N)\n",
    "        for i in np.arange(0, N - k):\n",
    "            autoCov += ((Xi[i + k]) - Xs) * (Xi[i] - Xs)\n",
    "        return (1 / (T)) * autoCov\n",
    "\n",
    "    gamma = []\n",
    "    for lag in range(0, h):\n",
    "        gamma.append(autocovariance(d_lst, len(d_lst), lag, mean_d))  # 0, 1, 2\n",
    "    V_d = (gamma[0] + 2 * sum(gamma[1:])) / T\n",
    "    DM_stat = V_d ** (-0.5) * mean_d\n",
    "    harvey_adj = ((T + 1 - 2 * h + h * (h - 1) / T) / T) ** (0.5)\n",
    "    DM_stat = harvey_adj * DM_stat\n",
    "    # Find p-value\n",
    "    p_value = 2 * t.cdf(-abs(DM_stat), df=T - 1)\n",
    "\n",
    "    # Construct named tuple for return\n",
    "    dm_return = collections.namedtuple('dm_return', 'DM p_value')\n",
    "    rt = dm_return(DM=DM_stat, p_value=p_value)\n",
    "    return rt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ETO</th>\n",
       "      <th>XGB</th>\n",
       "      <th>RF</th>\n",
       "      <th>CNN</th>\n",
       "      <th>GRU</th>\n",
       "      <th>MLP</th>\n",
       "      <th>Proposed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.243704</td>\n",
       "      <td>3.174300</td>\n",
       "      <td>3.792897</td>\n",
       "      <td>4.212972</td>\n",
       "      <td>3.943855</td>\n",
       "      <td>3.914148</td>\n",
       "      <td>3.290865</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.635035</td>\n",
       "      <td>2.466880</td>\n",
       "      <td>3.360410</td>\n",
       "      <td>3.118935</td>\n",
       "      <td>3.643656</td>\n",
       "      <td>3.079274</td>\n",
       "      <td>2.594209</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.537506</td>\n",
       "      <td>3.293312</td>\n",
       "      <td>3.393853</td>\n",
       "      <td>3.748376</td>\n",
       "      <td>2.853267</td>\n",
       "      <td>1.950341</td>\n",
       "      <td>3.319929</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.440828</td>\n",
       "      <td>3.457177</td>\n",
       "      <td>3.678498</td>\n",
       "      <td>3.336210</td>\n",
       "      <td>2.712084</td>\n",
       "      <td>2.057823</td>\n",
       "      <td>3.455395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.878740</td>\n",
       "      <td>3.976076</td>\n",
       "      <td>3.720689</td>\n",
       "      <td>3.634695</td>\n",
       "      <td>4.603443</td>\n",
       "      <td>4.256463</td>\n",
       "      <td>3.856467</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        ETO       XGB        RF       CNN       GRU       MLP  Proposed\n",
       "0  4.243704  3.174300  3.792897  4.212972  3.943855  3.914148  3.290865\n",
       "1  3.635035  2.466880  3.360410  3.118935  3.643656  3.079274  2.594209\n",
       "2  3.537506  3.293312  3.393853  3.748376  2.853267  1.950341  3.319929\n",
       "3  3.440828  3.457177  3.678498  3.336210  2.712084  2.057823  3.455395\n",
       "4  2.878740  3.976076  3.720689  3.634695  4.603443  4.256463  3.856467"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "data=pd.read_excel('../5 sub-models/寿县_28.xlsx')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [],
   "source": [
    "Observed=data['ETO']\n",
    "proposed=data['Proposed']\n",
    "lst1=data['XGB']\n",
    "lst2=data['RF']\n",
    "lst3=data['CNN']\n",
    "lst4=data['GRU']\n",
    "lst5=data['MLP']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "crit：MSE， dm_return(DM=-29.31587862385956, p_value=7.119324589400484e-174)\n"
     ]
    }
   ],
   "source": [
    "rt = dm_test(Observed, proposed, lst1, h=1, crit=\"MSE\")\n",
    "print(\"crit：MSE，\", rt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "crit：MSE， dm_return(DM=-7.705535247761928, p_value=1.5740255056879218e-14)\n"
     ]
    }
   ],
   "source": [
    "rt = dm_test(Observed, proposed, lst2, h=1, crit=\"MSE\")\n",
    "print(\"crit：MSE，\", rt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "crit：MSE， dm_return(DM=-0.8505104910055774, p_value=0.3950839305343212)\n"
     ]
    }
   ],
   "source": [
    "rt = dm_test(Observed, proposed, lst3, h=1, crit=\"MSE\")\n",
    "print(\"crit：MSE，\", rt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "crit：MSE， dm_return(DM=-3.028465785844647, p_value=0.002471114082885568)\n"
     ]
    }
   ],
   "source": [
    "rt = dm_test(Observed, proposed, lst4, h=1, crit=\"MSE\")\n",
    "print(\"crit：MSE，\", rt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "crit：MSE， dm_return(DM=-1.308238236516777, p_value=0.19085542070656586)\n"
     ]
    }
   ],
   "source": [
    "rt = dm_test(Observed, proposed, lst5, h=1, crit=\"MSE\")\n",
    "print(\"crit：MSE，\", rt)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DM<0  模型1比模型2更好                      DM>0  模型2比模型1更好"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "cddbdc538c09b0b47371ae40580457bfbde3af2437c2544403f3f5585f13c5a0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}