[PYTHON] [Introduction to infectious disease model] Looking at the logarithmic graph. .. .. One week from "the second wave"!

It's been another week since the last article, but I'd like to summarize the trends since then. This time, I would like to see the future situation using MACD, which is used for stock price forecasting. First, looking at the current situation on an ordinary logarithmic graph, Tokyo and the whole country are as follows. Tokyo has already exceeded the maximum value of the first wave and is still increasing exponentially. However, although it showed a slight downward trend from Sunday to Wednesday last week, it has no effect on the situation. Also, today is Sunday with 206 people, but the number of infections is almost the same. On the other hand, the whole country is also increasing exponentially. The whole country also showed a decrease from Sunday to Wednesday. The number of infections nationwide is almost the same as 371 today.

Apply MACD

First, it applies to the change of the first wave from the end of March to the end of April. The logic is briefly described as follows

MACD = 12-day EMA-26-day EMA
Signal = MACD 9-day EMA
Histogram = Signal-MACD

Here, EMA (Exponential Moving Average) is calculated by the following recurrence formula.

{S_{{t}}=\alpha \times Y_{{t-1}}+(1-\alpha )\times S_{{t-1}}
}

Here, $ Y_ {{t-1}} $ is the observed value at t-1 and $ S_ {{t-1}} $ is the previous EMA. $ \ Alpha = {2 \ over {N + 1}} $, N = 9 is a value of about 0.2, and $ S_ {{t}} $ is evaluated larger and averaged closer to the current observation point. Incorporated in. Therefore, the movement during the first wave can be calculated as follows. As shown below, the downward trend in Tokyo can be clearly read after April 17, but there is a lot of noise. The whole country is as follows, and almost the same as in Tokyo, a declining trend can be seen after April 17th. Actually, if you look at the above, the daily fluctuations are still large, so [Similar to the article on Systre, introduce Decompose here and use the trend curve for the same processing](https://qiita.com/MuAuan/items/ b08616a841be25d29817) is shown. Tokyo is as follows, and it can be seen that it has clearly turned to a declining trend after peaking on April 15. The whole country is as follows. This also peaked on April 16 and has turned to a declining trend. Therefore, it seems that this analysis method can be used to see the peak number of infections even in the time series of the number of infections in corona.

Apply to the second wave

The results in Tokyo are as follows. First, an analysis for no treatment; an exponential increase is also apparent in this analysis, with the number of infections likely to reach the level of 500 in 10 days to 2 weeks. Decompose and analyze the trend curve; there is no sign of further decline in this graph, and this trend is expected to continue for the time being. In other words, the MACD curve and signal curve are convex downward, suggesting that the difference will increase for the time being. In stock prices, the place where this difference turns negative can be regarded as a signal of decline, but this figure does not show that. The analysis of the whole country is as follows, and the curve is almost the same as that of Tokyo. Similarly, the analysis for the trend curve after decompose is as follows.

A little background analysis

The nationwide trend is similar to that of Tokyo. Of course, since the numbers in the Greater Tokyo area are large, it seems natural that the same tendency will occur, but the start of the second wave this time is June 22 from the slope of the blue plot and the rise of the signal-MACD bar graph. It is around. So, let's think about what causes this trigger. There are the following two trigger candidates. ① The state of emergency is canceled on May 25th. ② Tokyo alert cancellation is June 11th Comparing these, it is thought that it takes about 10 days from infection to detection by inspection, so it seems reasonable to think that the cancellation of the Tokyo alert is the trigger for the increase. To put it the other way around, it can be said that the Tokyo Alert had a certain effect.

Being concerned about the increase in the number of infections in the future

First of all, if you look at the red plot of the number of infections in Tokyo shown in the first graph, it has already increased to about 1200, and it seems that it will increase by the same index in the future. In other words, 3000 beds are likely to be filled in less than 10 days. It is reported that the usual length of hospital stay is about 10 days, so 3000 beds in less than 10 days is a little scary number. Next is the increase in critically ill patients. The increase in critically ill patients has been delayed and appears to be decreasing as the first wave of serious injuries are now recovering. However, this has begun to occur 7 to 10 days later in the first wave, and there is concern that it will occur in the future. And now, I hear that 70-80% of young people are, so the effect is suppressing the aggravation. However, as the number of infections increases, the absolute number of elderly people also increases, and the increase in the absolute number of aggravations cannot be prevented. And finally, 40 days after admission, when medical care begins to tighten, deaths begin to increase.

At this rate, there is no sign that the spread of infection will stop not only in Tokyo but throughout the country. We think it is necessary to take measures to prevent infection to high-risk individuals. As much as possible, it is assumed that the expansion cannot be suppressed unless it is implemented with the same strong driving force as the Tokyo Alert.

Summary

・ I applied MACD to predict the number of infections. ・ Exponential continuous increase was suggested ・ Anyway, the initial action (although it seems to be delayed) is important, so Tokyo Alert-class infection containment measures should be announced.

・ It is important what kind of countermeasure set the above alert should be while turning the economy, so I would like you to propose and implement it as soon as possible.

Bonus (code for the time being)

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt
from pandas_datareader import data
import statsmodels.api as sm
from statsmodels.tsa.seasonal import STL

#Read CSV data with pandas. C:\Users\user\simulation\COVID-19\csse_covid_19_data\japan\test
data = pd.read_csv('data/covid19/test_confirmed.csv',encoding="cp932")
data_r = pd.read_csv('data/covid19/test_recovered.csv',encoding="cp932")
data_d = pd.read_csv('data/covid19/test_deaths.csv',encoding="cp932")
#data = pd.read_csv('data/test_confirmed_.csv',encoding="cp932")
#data_r = pd.read_csv('data/test_recovered_.csv',encoding="cp932")
#data_d = pd.read_csv('data/test_deaths_.csv',encoding="cp932")

confirmed = [0] * (len(data.columns) - 1)
day_confirmed = [0] * (len(data.columns) - 1)
confirmed_r = [0] * (len(data_r.columns) - 1)
day_confirmed_r = [0] * (len(data.columns) - 1)
confirmed_d = [0] * (len(data_d.columns) - 1)
diff_confirmed = [0] * (len(data.columns) - 1)
days_from_1_Jun_20 = np.arange(0, len(data.columns) - 1, 1)
beta_ = [0] * (len(data_r.columns) - 1)
gamma_ = [0] * (len(data_d.columns) - 1)
daystamp = "531"

#city,city0 = "Tokyo","tokyo"
city,city0 = "total","total_japan"

skd=1
#Process the data
t_cases = 0
t_recover = 0
t_deaths = 0
for i in range(0, len(data_r), 1):
    if (data_r.iloc[i][0] == city): #for country/region
        print(str(data_r.iloc[i][0]))
        for day in range(1, len(data.columns), 1):            
            confirmed_r[day - 1] += data_r.iloc[i][day]
            if day < 1+skd:
                day_confirmed_r[day-1] += data_r.iloc[i][day]
            else:
                day_confirmed_r[day-1] += (data_r.iloc[i][day] - data_r.iloc[i][day-skd])/(skd)
        t_recover += data_r.iloc[i][day]        
for i in range(0, len(data_d), 1):
    if (data_d.iloc[i][0] == city): #for country/region
        print(str(data_d.iloc[i][0]) )
        for day in range(1, len(data.columns), 1):
            confirmed_d[day - 1] += float(data_d.iloc[i][day]) #fro drawings
        t_deaths += float(data_d.iloc[i][day])        
for i in range(0, len(data), 1):
    if (data.iloc[i][0] == city): #for country/region
        print(str(data.iloc[i][0]))
        for day in range(1, len(data.columns), 1):
            confirmed[day - 1] += data.iloc[i][day] -  confirmed_r[day - 1] -confirmed_d[day-1]
            if day == 1:
                day_confirmed[day-1] += data.iloc[i][day]
            else:
                day_confirmed[day-1] += data.iloc[i][day] - data.iloc[i][day-1]

def EMA1(x, n):
    a= 2/(n+1)
    return pd.Series(x).ewm(alpha=a).mean()                
                
day_confirmed[0]=0                
df = pd.DataFrame()                
date = pd.date_range("20200531", periods=len(day_confirmed))
df = pd.DataFrame(df,index = date)
df['Close'] = day_confirmed
df.to_csv('data/day_comfirmed_new_{}.csv'.format(city0))
date_df=df['Close'].index.tolist() #This is the point
print(date_df[0:30])
series = df['Close'].values.tolist()
stock0 = city0
stock = stock0
start = dt.date(2020,6,1)
end = dt.date(2020,7,12)

bunseki = "trend" #series" #cycle" #trend
cycle, trend = sm.tsa.filters.hpfilter(series, 144)
series2 = trend

y12 = EMA1(series2, 12)
y26 = EMA1(series2, 26)
MACD = y12 -y26
signal = EMA1(MACD, 9)
hist_=MACD-signal

ind3=date_df[:]
print(len(series),len(ind3))

fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 8, 4*2),dpi=200)
ax1.bar(ind3,series,label="series")
ax1.plot(ind3,series2, "o-", color="blue",label="series2")
ax1.plot(ind3,y12, ".-", color="red",label="y12")
ax1.plot(ind3,y26, ".-", color="green",label="y26")
ax2.plot(ind3,MACD,label="MACD")
ax2.plot(ind3,signal,label="signal")
ax2.bar(ind3,hist_)
ax1.legend()
ax2.legend()
ax1.set_ylim(10,1000)
#ax2.set_ylim(10,1000)
ax1.set_yscale('log')
#ax2.set_yscale('log')
ax1.grid()
ax2.grid()
plt.savefig("./fig/{}/ema_decompose_%5K%25D_{}_{}new{}.png ".format(stock0,stock,bunseki,start))
plt.pause(1)
plt.close()

df['Close']=series  #series" #cycle" #trend
df['series2']=series2
df['y12'] = EMA1(df['Close'], 12)
df['y26'] =  EMA1(df['Close'], 26)
df['MACD'] = df['y12'] -df['y26']
df['signal'] = EMA1(df['MACD'], 9)
df['hist_']=df['MACD']-df['signal']
date_df=df['Close'].index.tolist()
print(df[0:30])

fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 8, 4*2),dpi=200)
ax1.bar(ind3,series, label="series")
ax1.plot(df['series2'],"o-", color="blue",label="series2")
ax1.plot(df['y12'],".-", color="red",label="y12")
ax1.plot(df['y26'],".-", color="green",label="y26")
ax2.plot(df['MACD'],label="MACD")
ax2.plot(df['signal'],label="signal")
ax2.bar(date_df,df['hist_'])
ax1.legend()
ax2.legend()
ax1.set_ylim(10,1000)
ax1.set_yscale('log')
ax1.grid()
ax2.grid()
plt.savefig("./fig/{}/ema_df_decompose_%5K%25D_{}_{}new{}.png ".format(stock0,stock,bunseki,start))
plt.pause(1)
plt.close()