{
 "cells": [
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T15:31:07.837207Z",
     "start_time": "2025-05-23T15:31:05.479521Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.tsa.ar_model import AutoReg\n",
    "from statsmodels.tsa.stattools import adfuller\n",
    "from statsmodels.graphics.tsaplots import plot_pacf\n",
    "from statsmodels.tsa.ar_model import ar_select_order"
   ],
   "id": "3d3cd073f42a450a",
   "outputs": [],
   "execution_count": 1
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T15:31:07.881789Z",
     "start_time": "2025-05-23T15:31:07.866734Z"
    }
   },
   "cell_type": "code",
   "source": [
    "#import data\n",
    "main_folder = Path.cwd().parent.parent\n",
    "\n",
    "data_folder = main_folder / 'data/Norbert'\n",
    "plots_folder = main_folder / 'Plots/Norbert'\n",
    "scripts_folder = main_folder / 'scripts/Norbert'\n",
    "\n",
    "df = pd.read_csv(data_folder / 'deaths_malaysia.csv')"
   ],
   "id": "initial_id",
   "outputs": [],
   "execution_count": 2
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T15:31:09.768313Z",
     "start_time": "2025-05-23T15:31:09.738628Z"
    }
   },
   "cell_type": "code",
   "source": [
    "df.drop(columns = ['deaths_new', 'deaths_bid', 'deaths_new_dod', 'deaths_bid_dod', 'deaths_pvax', 'deaths_boost', 'deaths_tat'], inplace = True)\n",
    "df.columns = ['date', 'unvax', 'fvax']\n",
    "df.set_index('date', inplace = True)\n",
    "df.index = pd.to_datetime(df.index)\n",
    "df = df.asfreq('D')\n",
    "print(df)\n"
   ],
   "id": "3810b31b555cc0",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            unvax  fvax\n",
      "date                   \n",
      "2020-03-17      2     0\n",
      "2020-03-18      0     0\n",
      "2020-03-19      0     0\n",
      "2020-03-20      1     0\n",
      "2020-03-21      4     0\n",
      "...           ...   ...\n",
      "2025-04-29      0     0\n",
      "2025-04-30      0     0\n",
      "2025-05-01      0     0\n",
      "2025-05-02      0     0\n",
      "2025-05-03      0     0\n",
      "\n",
      "[1874 rows x 2 columns]\n"
     ]
    }
   ],
   "execution_count": 3
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T14:44:36.849829Z",
     "start_time": "2025-05-23T14:44:36.600068Z"
    }
   },
   "cell_type": "code",
   "source": [
    "#check data stationarity\n",
    "adf_result = adfuller(df['unvax'].dropna())\n",
    "print('ADF Statistic for unvax:', adf_result[0])  # More negative = stronger evidence for stationarity\n",
    "print('p-value for unvax:', adf_result[1])\n",
    "\n",
    "adf_result = adfuller(df['fvax'].dropna())\n",
    "print('ADF Statistic for fvax:', adf_result[0])  # More negative = stronger evidence for stationarity\n",
    "print('p-value for fvax:', adf_result[1])\n",
    "\n",
    "pcaf_unvax = plot_pacf(df['unvax'], lags = 30)\n",
    "pcaf_fvax = plot_pacf(df['fvax'], lags = 30)\n",
    "\n",
    "\n",
    "model_unvax = AutoReg(df['unvax'], lags = 4, old_names=False)\n",
    "model_unvax_fitted = model_unvax.fit()\n",
    "print(model_unvax_fitted.summary())\n",
    "\n",
    "model_fvax = AutoReg(df['fvax'], lags = 4, old_names=False)\n",
    "model_fvax_fitted = model_fvax.fit()\n",
    "print(model_fvax_fitted.summary())\n"
   ],
   "id": "4fd784f14eac254d",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ADF Statistic for unvax: -3.0963655667430716\n",
      "p-value for unvax: 0.026823324581619963\n",
      "ADF Statistic for fvax: -3.5181435549832507\n",
      "p-value for fvax: 0.0075331737446544915\n",
      "                            AutoReg Model Results                             \n",
      "==============================================================================\n",
      "Dep. Variable:                  unvax   No. Observations:                 1874\n",
      "Model:                     AutoReg(4)   Log Likelihood               -5649.105\n",
      "Method:               Conditional MLE   S.D. of innovations              4.963\n",
      "Date:                Fri, 23 May 2025   AIC                          11310.210\n",
      "Time:                        16:44:36   BIC                          11343.412\n",
      "Sample:                    03-21-2020   HQIC                         11322.443\n",
      "                         - 05-03-2025                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.0594      0.122      0.488      0.625      -0.179       0.298\n",
      "unvax.L1       0.4122      0.023     18.146      0.000       0.368       0.457\n",
      "unvax.L2       0.2649      0.024     10.836      0.000       0.217       0.313\n",
      "unvax.L3       0.1311      0.024      5.363      0.000       0.083       0.179\n",
      "unvax.L4       0.1867      0.023      8.219      0.000       0.142       0.231\n",
      "                                    Roots                                    \n",
      "=============================================================================\n",
      "                  Real          Imaginary           Modulus         Frequency\n",
      "-----------------------------------------------------------------------------\n",
      "AR.1            1.0024           -0.0000j            1.0024           -0.0000\n",
      "AR.2           -1.7065           -0.0000j            1.7065           -0.5000\n",
      "AR.3            0.0010           -1.7694j            1.7694           -0.2499\n",
      "AR.4            0.0010           +1.7694j            1.7694            0.2499\n",
      "-----------------------------------------------------------------------------\n",
      "                            AutoReg Model Results                             \n",
      "==============================================================================\n",
      "Dep. Variable:                   fvax   No. Observations:                 1874\n",
      "Model:                     AutoReg(4)   Log Likelihood               -4271.492\n",
      "Method:               Conditional MLE   S.D. of innovations              2.376\n",
      "Date:                Fri, 23 May 2025   AIC                           8554.984\n",
      "Time:                        16:44:36   BIC                           8588.186\n",
      "Sample:                    03-21-2020   HQIC                          8567.217\n",
      "                         - 05-03-2025                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.0487      0.059      0.823      0.411      -0.067       0.165\n",
      "fvax.L1        0.2431      0.022     10.826      0.000       0.199       0.287\n",
      "fvax.L2        0.2631      0.022     11.711      0.000       0.219       0.307\n",
      "fvax.L3        0.2422      0.022     10.784      0.000       0.198       0.286\n",
      "fvax.L4        0.2386      0.022     10.624      0.000       0.195       0.283\n",
      "                                    Roots                                    \n",
      "=============================================================================\n",
      "                  Real          Imaginary           Modulus         Frequency\n",
      "-----------------------------------------------------------------------------\n",
      "AR.1            1.0053           -0.0000j            1.0053           -0.0000\n",
      "AR.2           -1.6516           -0.0000j            1.6516           -0.5000\n",
      "AR.3           -0.1845           -1.5781j            1.5889           -0.2685\n",
      "AR.4           -0.1845           +1.5781j            1.5889            0.2685\n",
      "-----------------------------------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "[encoded data, 23140 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "[encoded data, 23752 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 262
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T14:16:39.029363Z",
     "start_time": "2025-05-23T14:16:38.981179Z"
    }
   },
   "cell_type": "code",
   "source": [
    "forecast_steps = 3652\n",
    "forecast_diff_unvax = model_unvax_fitted.predict(start = len(df['unvax']), end = len(df['unvax']) + forecast_steps-1, dynamic = False)\n",
    "forecast_diff_fvax = model_fvax_fitted.predict(start = len(df['fvax']), end = len(df['fvax']) + forecast_steps-1, dynamic = False)\n",
    "\n",
    "last_deaths_unvax = df['unvax'].iloc[-1]\n",
    "last_deaths_fvax = df['fvax'].iloc[-1]\n",
    "\n",
    "forecast_deaths_unvax = last_deaths_unvax + forecast_diff_unvax\n",
    "forecast_deaths_fvax = last_deaths_fvax + forecast_diff_fvax\n",
    "\n",
    "forecast_date = pd.date_range(start = df.index[-1] + pd.Timedelta(days =1), periods = forecast_steps, freq = 'D')\n",
    "\n",
    "forecast_df_unvax = pd.DataFrame({'data': forecast_date, 'deaths': forecast_deaths_unvax})\n",
    "forecast_df_fvax = pd.DataFrame({'data': forecast_date, 'deaths': forecast_deaths_fvax})"
   ],
   "id": "8f3e385cf6971019",
   "outputs": [],
   "execution_count": 232
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T14:41:40.511282Z",
     "start_time": "2025-05-23T14:41:40.251943Z"
    }
   },
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(df.index, df['unvax'],  label = 'Unvaccinated', color = '0.1')\n",
    "plt.plot(df.index, df['fvax'],label = 'Fully Vaccinated', color = '0.7')\n",
    "plt.plot(forecast_df_unvax.index, forecast_df_unvax['deaths'], color = '0.2')\n",
    "plt.plot(forecast_df_fvax.index, forecast_df_fvax['deaths'], color = '0.6')\n",
    "plt.axvline(x = df.index[-1], linestyle = '--', color = 'black')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Deaths Daily')\n",
    "plt.title('Daily Death Predictions across fully vaccinated and unvaccinated Malyasia citizens')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(plots_folder/'daily_deaths.jpg')\n",
    "\n"
   ],
   "id": "a4fdc27d4944de12",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ],
      "image/png": "[encoded data, 48572 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 259
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T14:16:41.537500Z",
     "start_time": "2025-05-23T14:16:41.531709Z"
    }
   },
   "cell_type": "code",
   "source": [
    "last_deaths_unvax = df['unvax'].iloc[-1]\n",
    "last_deaths_fvax = df['fvax'].iloc[-1]\n",
    "\n",
    "forecast_deaths_unvax = last_deaths_unvax + forecast_diff_unvax\n",
    "forecast_deaths_fvax = last_deaths_fvax + forecast_diff_fvax\n",
    "\n",
    "forecast_date = pd.date_range(start=df.index[-1] + pd.Timedelta(days=1), periods=forecast_steps, freq='D')\n",
    "\n",
    "forecast_df_unvax = pd.DataFrame({'data': forecast_date, 'deaths': forecast_deaths_unvax})\n",
    "forecast_df_fvax = pd.DataFrame({'data': forecast_date, 'deaths': forecast_deaths_fvax})"
   ],
   "id": "c6f4a17c6e28e827",
   "outputs": [],
   "execution_count": 235
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-05-23T14:40:57.260731Z",
     "start_time": "2025-05-23T14:40:57.046509Z"
    }
   },
   "cell_type": "code",
   "source": [
    "cum_df = df.cumsum()\n",
    "last_cum_deaths_unvax = cum_df['unvax'].iloc[-1]\n",
    "last_cum_deaths_fvax = cum_df['fvax'].iloc[-1]\n",
    "\n",
    "forecast_deaths_cum_unvax = last_cum_deaths_unvax + forecast_diff_unvax.cumsum()\n",
    "forecast_deaths_cum_fvax = last_cum_deaths_fvax + forecast_diff_fvax.cumsum()\n",
    "forecast_date = pd.date_range(start = df.index[-1] + pd.Timedelta(days =1), periods = forecast_steps, freq = 'D')\n",
    "\n",
    "forecast_df_deaths_cum_unvax = pd.DataFrame({'data': forecast_date, 'deaths': forecast_deaths_cum_unvax})\n",
    "forecast_df_deaths_cum_fvax = pd.DataFrame({'data': forecast_date, 'deaths': forecast_deaths_cum_fvax})\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(cum_df.index, cum_df['unvax'], label = 'Unvaccinated', color = '0.1')\n",
    "plt.plot(cum_df.index, cum_df['fvax'], label = 'Fully Vaccinated', color = '0.7')\n",
    "plt.plot(forecast_df_unvax.index, forecast_df_deaths_cum_unvax['deaths'], color = '0.2')\n",
    "plt.plot(forecast_df_fvax.index, forecast_df_deaths_cum_fvax['deaths'], color = '0.6')\n",
    "plt.axvline(x = cum_df.index[-1], linestyle = '--', color = 'black', label = 'Prediction start')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Cumulative Deaths')\n",
    "plt.title('Cumulative COVID-19 Death Forecast: Fully Vaccinated vs. Unvaccinated in Malaysia')\n",
    "plt.tight_layout()\n",
    "plt.legend()\n",
    "plt.savefig(plots_folder/'cumulative_deaths.jpg')\n"
   ],
   "id": "d73f1691e16e32b8",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ],
      "image/png": "[encoded data, 65180 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 256
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": "",
   "id": "b95c32d455df800a"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
