{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Quick study to investigate oscillations in reported infections in Germany. Here is the plot of the data in question: " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/markdown": [ "[Execute this notebook with Binder](https://mybinder.org/v2/gh/oscovida/binder/master?filepath=ipynb/2020-05-10-notebook-weekly-fluctuations-in-data-from-germany.ipynb)" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import coronavirus\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%config InlineBackend.figure_formats = ['svg']\n", "\n", "coronavirus.display_binder_link(\"2020-05-10-notebook-weekly-fluctuations-in-data-from-germany.ipynb\")" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# get data\n", "cases, deaths, country_label = coronavirus.get_country_data(\"Germany\")\n", "\n", "# plot daily changes\n", "fig, ax = plt.subplots(figsize=(8, 4))\n", "coronavirus.plot_daily_change(ax, cases, 'C1')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The working assumption is that during the weekend fewer numbers are captured or reported. The analysis below seems to confirm this.\n", "\n", "We compute a discrete Fourier transform of the data, and expect a peak at a frequency corresponding to a period of 7 days." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data selection\n", "\n", "We start with data from 1st March as numbers before were small. It is convenient to take a number of days that can be divided by seven (for alignment of the freuency axis in Fourier space, so we choose 63 days from 1st of March):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "63" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data = cases['2020-03-01':'2020-05-03']\n", "# compute daily change\n", "diff = data.diff().dropna()\n", "# plot data points (corresponding to bars in figure above:)\n", "fig, ax = plt.subplots()\n", "ax.plot(diff.index, diff, '-C1', \n", " label='daily new cases Germany')\n", "\n", "fig.autofmt_xdate() # avoid x-labels overlap\n", "# How many data points (=days) have we got?\n", "diff.size" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "63" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "diff2 = diff.resample(\"24h\").asfreq() # ensure we have one data point every day\n", "diff2.size" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute the frequency spectrum" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "# compute power density spectrum \n", "change_F = abs(np.fft.fft(diff2))**2\n", "# determine appropriate frequencies\n", "n = change_F.size\n", "freq = np.fft.fftfreq(n, d=1)\n", "# We skip values at indices 0, 1 and 2: these are large because we have a finite \n", "# sequence and not substracted the mean from the data set\n", "# We also only plot the the first n/2 frequencies as for high n, we get negative\n", "# frequencies with the same data content as the positive ones.\n", "ax.plot(freq[3:n//2], change_F[3:n//2], 'o-C3')\n", "ax.set_xlabel('frequency [cycles per day]');\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A signal with oscillations on a weekly basis would correspond to a frequency of 1/7 as frequency is measured in `per day`. We thus expect the peak above to be at 1/7 $\\approx 0.1428$. \n", "\n", "We can show this more easily by changing the frequency scale from cycles per day to cycles per week:\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(freq[3:n//2] * 7, change_F[3:n//2], 'o-C3')\n", "ax.set_xlabel('frequency [cycles per week]');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In other words: there as a strong component of the data with a frequency corresponding to one week.\n", "\n", "\n", "This is the end of the notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Fourier transform basics\n", "\n", "A little playground to explore properties of discrete Fourier transforms." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "time = np.linspace(0, 4, 1000)\n", "signal_frequency = 3 # choose this freely\n", "signal = np.sin(time * 2 * np.pi * signal_frequency)\n", "fourier = np.abs(np.fft.fft(signal))\n", "\n", "# compute frequencies in fourier spectrum\n", "n = signal.size\n", "timestep = time[1] - time[0]\n", "freqs = np.fft.fftfreq(n, d=timestep)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(time, signal, 'oC9', label=f'signal, frequency={signal_frequency}')\n", "ax.set_xlabel('time')\n", "ax.legend()\n", "fig, ax = plt.subplots()\n", "ax.plot(freqs[0:n//2][:20], fourier[0:n//2][0:20], 'o-C8', label=\"Fourier transform\")\n", "ax.legend()\n", "ax.set_xlabel('frequency');\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/markdown": [ "[Execute this notebook with Binder](https://mybinder.org/v2/gh/oscovida/binder/master?filepath=ipynb/2020-05-10-notebook-weekly-fluctuations-in-data-from-germany.ipynb)" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "coronavirus.display_binder_link(\"2020-05-10-notebook-weekly-fluctuations-in-data-from-germany.ipynb\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }