Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data were used is available in [Database source](https://github.com/owid/covid-19-data/tree/master/public/data)\n",
    "\n",
    "#### Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy import asarray as ar,exp\n",
    "import math\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats as stats\n",
    "import scipy.linalg as linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading Data\n",
    "\n",
    "Loading data fill NA with 0 and get all contries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = pd.read_csv('world.csv',sep=',')\n",
    "world.fillna(0,inplace=True)\n",
    "countries = world['location'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose the number of days to create a gaussian model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "begin = 1\n",
    "end   = 150"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variable for configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphicDays = 200\n",
    "factor = 40\n",
    "significance = 0.05\n",
    "directory = 'calculated/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Referência para o p-value a partir do F](https://www.programcreek.com/python/example/113313/scipy.stats.f.cdf)\n",
    "\n",
    "### Returns the p-value using the computed F-statistic\n",
    "\n",
    "#### Returns\n",
    "\n",
    "p : float\n",
    "    The computed p-value given the found F-statistic and the group and residual degrees of freedom.\n",
    "\n",
    "#### Notes\n",
    "\n",
    "The code `cdf` method from Scipy's code `stats.f` class is used to find the p-value.\n",
    "\n",
    "#### References\n",
    "\n",
    "Rencher, A. (n.d.). Methods of Multivariate Analysis (2nd ed.). Brigham Young University: John Wiley & Sons, Inc.\n",
    "\n",
    "#### Códigos\n",
    "`\n",
    "p = 1 - f.cdf(self.f_statistic,\n",
    "      self.group_degrees_of_freedom,\n",
    "      self.residual_degrees_of_freedom)\n",
    "`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Used model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_fun(x,a,b,c):\n",
    "    return  a * np.exp(-(x - b) ** 2 / (2 * c ** 2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chow's test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chow_test(y1,y2):\n",
    "    x1 = ar(list(range(len(y1))))\n",
    "    x2 = ar(list(range(len(y2))))\n",
    "    x3 = ar(list(range(len(y1)+len(y2))))\n",
    "    y3 = np.append(y1,y2)\n",
    "    \n",
    "    x1 = sm.add_constant(x1)\n",
    "    x2 = sm.add_constant(x2)\n",
    "    x3 = sm.add_constant(x3)\n",
    "    \n",
    "    model1 = sm.OLS(y1,x1).fit()\n",
    "    model2 = sm.OLS(y2,x2).fit()\n",
    "    model3 = sm.OLS(y3,x3).fit()\n",
    "    \n",
    "    RSSall = model3.ssr\n",
    "    RSSbefore = model1.ssr\n",
    "    RSSafter = model2.ssr\n",
    "    \n",
    "    N1 = len(x1)\n",
    "    N2 = len(x2)\n",
    "    chow = ((RSSall-(RSSbefore+RSSafter))/2.0)/((RSSbefore+RSSafter)/(N1+N2-4.0))\n",
    "    \n",
    "    df1 = 2.0\n",
    "    df2 = len(x1) + len(x2) - 4.0\n",
    " \n",
    "    p_val = 1 - stats.f.cdf(chow, df1, df2)\n",
    "    return p_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/raphael/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py:787: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Afghanistan 3.4467764438101263e-10\n",
      "Africa 6.593652577808076e-10\n",
      "Albania 5.869709006343797e-07\n",
      "Algeria 5.181167722702819e-06\n",
      "Andorra fora\n",
      "Angola 0.0001557813624260258\n",
      "Antigua and Barbuda 0.0007516421120277423\n",
      "Argentina 7.725660463092637e-09\n",
      "Armenia 2.7061698891566688e-08\n",
      "Asia 1.2384507307578361e-09\n",
      "Australia 6.646083749602794e-05\n",
      "Austria fora\n",
      "Azerbaijan 1.4463192508568895e-05\n",
      "Bahamas fora\n",
      "Bahrain 8.194455522433007e-05\n",
      "Bangladesh fora\n",
      "Barbados fora\n",
      "Belarus 8.285897956459051e-12\n",
      "Belgium 1.1583874822269552e-06\n",
      "Belize fora\n",
      "Benin 2.3099599147610657e-10\n",
      "Bolivia 1.2925555800900236e-07\n",
      "Bosnia and Herzegovina 3.925454948169911e-09\n",
      "Botswana 0.28999723824619755\n",
      "Brazil 7.24878343073568e-11\n",
      "Brunei fora\n",
      "Bulgaria 5.310359070408066e-05\n",
      "Burkina Faso fora\n",
      "Burundi 0.526690736482871\n",
      "Cambodia fora\n",
      "Cameroon 6.387934202427651e-07\n",
      "Canada 0.000505858707614583\n",
      "Cape Verde 2.021047647265454e-05\n",
      "Central African Republic 4.01251998509658e-07\n",
      "Chad fora\n",
      "Chile 1.1303409175848528e-08\n",
      "China 1.2131635873850863e-09\n",
      "Colombia 3.3794219321620463e-10\n",
      "Comoros fora\n",
      "Congo 0.003532996667467553\n",
      "Costa Rica 8.75422153362579e-11\n",
      "Cote d'Ivoire fora\n",
      "Croatia 5.734512862489635e-08\n",
      "Cuba fora\n",
      "Cyprus fora\n",
      "Czechia fora\n",
      "Democratic Republic of Congo 0.00018665781838101702\n",
      "Denmark 0.06353096964705787\n",
      "Djibouti 2.2945219060630428e-11\n",
      "Dominica fora\n",
      "Dominican Republic 1.465309327763044e-09\n",
      "Ecuador 0.028699122486153125\n",
      "Egypt 9.807388325024669e-12\n",
      "El Salvador 7.115575689330788e-11\n",
      "Equatorial Guinea 0.06422587901634795\n",
      "Eritrea fora\n",
      "Estonia fora\n",
      "Eswatini 0.0043501254665555406\n",
      "Ethiopia 2.917123618952146e-12\n",
      "Europe 0.00022717675388105927\n",
      "European Union 0.0001920133480672842\n",
      "Finland fora\n",
      "France 3.726582275562577e-05\n",
      "Gabon 7.036049807122385e-07\n",
      "Gambia 6.275051110742367e-11\n",
      "Georgia 2.7675810445166037e-05\n",
      "Germany 0.0008728309793441546\n",
      "Ghana 3.5833450318998624e-05\n",
      "Greece fora\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:1362: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return self.params / self.bse\n",
      "/home/raphael/.local/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:877: RuntimeWarning: invalid value encountered in greater\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/home/raphael/.local/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:877: RuntimeWarning: invalid value encountered in less\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/home/raphael/.local/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:1831: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= self.a)\n",
      "/usr/lib/python3/dist-packages/ipykernel_launcher.py:21: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Grenada nan\n",
      "Guatemala 1.1225406668639962e-09\n",
      "Guinea fora\n",
      "Guinea-Bissau fora\n",
      "Guyana 0.00011056736436210619\n",
      "Haiti 0.0014819075723233866\n",
      "Honduras 0.011509380624530646\n",
      "Hungary fora\n",
      "Iceland 2.8050534028008257e-08\n",
      "India 2.0872626007223607e-05\n",
      "Indonesia 1.848457560546312e-05\n",
      "International fora\n",
      "Iran 1.8752947905269558e-08\n",
      "Iraq 1.3917808661651145e-05\n",
      "Ireland 1.764586180709739e-05\n",
      "Israel 2.7813452467245294e-11\n",
      "Italy 1.2702530587267272e-07\n",
      "Jamaica fora\n",
      "Japan 1.548807526478338e-08\n",
      "Jordan 0.09711711638074558\n",
      "Kazakhstan 3.975581505590302e-09\n",
      "Kenya 2.905430778103939e-07\n",
      "Kosovo fora\n",
      "Kuwait 1.5123550568525951e-06\n",
      "Kyrgyzstan 0.38157864482484805\n",
      "Latvia fora\n",
      "Lebanon 0.2850499129854625\n",
      "Lesotho fora\n",
      "Liberia 0.047084651250427115\n",
      "Libya 8.625899708369656e-10\n",
      "Liechtenstein fora\n",
      "Lithuania fora\n",
      "Luxembourg fora\n",
      "Madagascar 0.00013946662486873408\n",
      "Malawi 0.01470760285593435\n",
      "Malaysia fora\n",
      "Maldives 1.5657296883538816e-05\n",
      "Mali fora\n",
      "Malta fora\n",
      "Mauritania 1.8267141560344822e-08\n",
      "Mauritius fora\n",
      "Mexico 4.6170039998398453e-10\n",
      "Moldova 0.0036021935594980972\n",
      "Monaco fora\n",
      "Mongolia fora\n",
      "Montenegro 1.8713901646649987e-06\n",
      "Morocco 1.7362648466358393e-11\n",
      "Mozambique 0.022804655999746848\n",
      "Myanmar 6.429852640343906e-07\n",
      "Namibia 1.393836790299273e-10\n",
      "Nepal 8.306572738805802e-08\n",
      "Netherlands fora\n",
      "New Zealand fora\n",
      "Nicaragua fora\n",
      "Niger fora\n",
      "Nigeria 3.7125760978246094e-08\n",
      "North America 0.0014968582055386599\n",
      "North Macedonia 0.18414516431361322\n",
      "Norway 8.533591552658923e-10\n",
      "Oceania 1.1656933136605397e-10\n",
      "Oman 5.602573438109113e-13\n",
      "Pakistan 4.0498582073770425e-08\n",
      "Palestine 3.322542488669748e-08\n",
      "Panama 4.041381461567527e-07\n",
      "Papua New Guinea 2.8304061469531196e-07\n",
      "Paraguay fora\n",
      "Peru 3.947396381182711e-05\n",
      "Philippines 5.1972135062937013e-11\n",
      "Poland fora\n",
      "Portugal fora\n",
      "Qatar 5.6601053030652154e-12\n",
      "Romania 5.3077335373221755e-09\n",
      "Russia 9.112223201910812e-05\n",
      "Rwanda 1.0500549697140061e-06\n",
      "Saint Kitts and Nevis fora\n",
      "San Marino fora\n",
      "Sao Tome and Principe fora\n",
      "Saudi Arabia 1.6478031473608535e-12\n",
      "Senegal fora\n",
      "Serbia 2.3502120014990236e-08\n",
      "Seychelles fora\n",
      "Sierra Leone fora\n",
      "Singapore fora\n",
      "Slovakia 9.431467454832715e-05\n",
      "Slovenia fora\n",
      "Somalia 7.092114971435274e-10\n",
      "South Africa 1.1160635069529558e-07\n",
      "South America 9.043864816263065e-09\n",
      "South Korea fora\n",
      "South Sudan fora\n",
      "Spain 1.292711255199255e-08\n",
      "Sri Lanka fora\n",
      "Sudan fora\n",
      "Suriname 1.0069845767118514e-08\n",
      "Sweden 6.033871303943912e-12\n",
      "Switzerland 7.79078177403607e-07\n",
      "Syria 0.002371213059490025\n",
      "Taiwan fora\n",
      "Tajikistan fora\n",
      "Tanzania fora\n",
      "Thailand 3.1682300626354667e-12\n",
      "Togo 8.916609779403094e-10\n",
      "Trinidad and Tobago fora\n",
      "Tunisia fora\n",
      "Turkey 1.3360954053433344e-09\n",
      "Uganda fora\n",
      "Ukraine 1.2417359332152634e-10\n",
      "United Arab Emirates 8.998067776746808e-06\n",
      "United Kingdom 7.336825143700021e-06\n",
      "United States 0.00490030696919906\n",
      "Uruguay 0.00026843798570553315\n",
      "Uzbekistan 1.5520332857793253e-07\n",
      "Venezuela 0.40065026102156653\n",
      "Vietnam fora\n",
      "World 1.2016948122430931e-09\n",
      "Yemen 0.017376520773387902\n",
      "Zambia 1.10262973275752e-06\n",
      "Zimbabwe 0.016378975694786357\n",
      "END 100 24 63 124 187 25 212\n"
     ]
    }
   ],
   "source": [
    "significant = 0\n",
    "nonSignificant = 0\n",
    "out = 0\n",
    "for country in countries:\n",
    "    significantDate = None\n",
    "    dataCountry = world[world['location'] == country]\n",
    "    density = dataCountry['population_density'].values[0]\n",
    "    if(density <= 0):\n",
    "        density = 25\n",
    "\n",
    "    means = dataCountry['new_cases'].rolling(7).mean()\n",
    "    means[:6] = dataCountry['new_cases'].values[:6]\n",
    "    dateX = dataCountry['date']\n",
    "    x = ar(range(len(means)))\n",
    "    y = ar(means.values)\n",
    "    \n",
    "\n",
    "    coefs = []\n",
    "    dates = []\n",
    "    endding = end\n",
    "    if(len(y) < endding):\n",
    "        endding = len(y)\n",
    "    getIn = False\n",
    "    for f in range(begin,endding):\n",
    "        if((y[f] > 2*density and y[f-1] > 2*density) or f >= 30):\n",
    "            trainX = x[:f]\n",
    "            trainY = y[:f]\n",
    "            if(trainY.max() > 3):\n",
    "                try:\n",
    "                    popt,pcov = curve_fit(model_fun,trainX,trainY,[1, 1, 100], maxfev=5000)\n",
    "                except:\n",
    "                    try:\n",
    "                        popt,pcov = curve_fit(model_fun,trainX,trainY,[10, 1, 1], maxfev=5000)\n",
    "                    except:\n",
    "                        continue\n",
    "            else:\n",
    "                continue\n",
    "            dates.append(dateX.values[f])\n",
    "            coefs.append(popt[0])\n",
    "            \n",
    "            if(significantDate == None and popt[0] >= factor*density and len(coefs) > 3 and f >= 14):\n",
    "                significantDate = dateX.values[f]\n",
    "                getIn = True\n",
    "                before = y[f-14:f]\n",
    "                after = y[f:f+14]\n",
    "\n",
    "                df = pd.DataFrame()\n",
    "                g1 = np.full(len(before),1)\n",
    "                g2 = np.full(len(after),2)\n",
    "                \n",
    "                d1 = list(zip(g1,before))\n",
    "                d2 = list(zip(g2, after))\n",
    "                \n",
    "                df = df.append(d1)\n",
    "                df = df.append(d2)\n",
    "                df.reset_index(drop=True, inplace=True)\n",
    "\n",
    "                df.columns = ['Group', 'Cases']\n",
    "                results = smf.ols('Cases ~ C(Group)', data=df).fit()\n",
    "\n",
    "                pchow = chow_test(before,after)\n",
    "                if(results.pvalues[1] <= significance and pchow <= significance):\n",
    "                    significant = significant + 1\n",
    "                else:\n",
    "                    nonSignificant = nonSignificant + 1\n",
    "\n",
    "    if(not getIn and len(coefs) > 0):\n",
    "        out = out + 1\n",
    "        \n",
    "    if(significantDate != None):\n",
    "        print(country)\n",
    "        dates = ar(dates)\n",
    "        fig = plt.figure(figsize=(16,9))\n",
    "        plt.plot(dates,coefs,'b+:',label='Gaussiana')\n",
    "        plt.yscale('log')\n",
    "        plt.title(country)\n",
    "        idx = list(range(0,len(dates),2))\n",
    "        plt.xticks(dates[idx], dates[idx], rotation='vertical')\n",
    "        plt.hlines(factor*density,0,len(coefs))\n",
    "        plt.legend()\n",
    "\n",
    "        try:\n",
    "            plt.savefig(directory+country+'_gaussian.png')\n",
    "        except:\n",
    "            print('ERRO:',country)\n",
    "            pass\n",
    "        fig.clf()\n",
    "        plt.close()\n",
    "    \n",
    "        fig = plt.figure(figsize=(16,9))\n",
    "        plt.plot(dateX[:graphicDays],y[:graphicDays],label='New Cases')\n",
    "\n",
    "        plt.vlines(significantDate,means[:graphicDays].min(),means[:graphicDays].max())\n",
    "        plt.title(country)\n",
    "        idx = list(range(0,len(dateX[:graphicDays]),10))\n",
    "        dateX = dateX.reset_index(drop=True)\n",
    "        plt.xticks(dateX[idx], dateX[idx], rotation='vertical')\n",
    "        plt.legend()\n",
    "\n",
    "        try:\n",
    "            plt.savefig(directory+country+'_curva.png')\n",
    "        except:\n",
    "            print('ERRO:',country)\n",
    "            pass\n",
    "\n",
    "        fig.clf()\n",
    "        plt.close()\n",
    "    #break\n",
    "    \n",
    "    \n",
    "print('END', significant, nonSignificant, out, significant+nonSignificant, significant+nonSignificant+out, len(countries)-significant-nonSignificant-out, len(countries))"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}